Very first thing first, to start the evaluation we have to have available the precise territory of those cities. To get them you may both use Google Earth Engine dataset known as FAO GAUL: World Administrative Unit Layers 2015 or the GADM web site. Consequently we must always find yourself with a bunch of polygons every resembling a Ukrainian area.
To create such a visualization you’ll must obtain the aforementioned boundaries and browse them utilizing geopandas library:
form = gpd.read_file('YOUR_FILE.shp')
form = form[(shape['NAME_1']=='Kiev') | (form['NAME_1']=='Kiev Metropolis') | (form['NAME_1']=='?') | (form['NAME_1']=='Kharkiv')|
(form['NAME_1']=='Odessa')]
form.plot(coloration='gray', edgecolor='black')plt.axis('off')
plt.textual content(35,48, 'Kharkov', fontsize=20)
plt.textual content(31,46, 'Odessa', fontsize=20)
plt.textual content(31,49, 'Kiev', fontsize=20)
plt.savefig('UKR_shape.png')
plt.present()
The second step for us goes to be acquisition of VIIRS information via GEE. For those who obtain the form of Ukrainian areas from the web, you’ll must wrap it right into a GEE geometry object. In any other case you’ve already received it prepared to make use of.
import json
import eejs = json.hundreds(form.to_json())
roi = ee.Geometry(ee.FeatureCollection(js).geometry())
Now let’s outline the timeline of our analysis. Conceptually, to know if the evening time gentle radiance after the start of the conflict was anomalous, we have to know the values earlier than. So we will probably be working with the entire timeframe out there: from 2012–01–01 to 2024–04–01. Knowledge earlier than 2022–02–01 will probably be thought of as «a norm» and every thing after will probably be subtracted from this norm, therefore, representing a deviation (anomaly).
startDate = pd.to_datetime('2012-01-01')
endDate = pd.to_datetime('2024-04-01')
information = ee.ImageCollection("NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG")
.filterBounds(roi)
.filterDate(begin = startDate, finish=endDate)
Our ultimate consequence will embrace a map and anomaly plot. To carry out this visualization we have to accumulate month-to-month evening time lights radiance maps between 2022–02–01 and 2024–04–01 and common month-to-month evening time lights radiance (within the type of time collection) for every area. One of the simplest ways to do this is to iterate via an inventory of GEE photos and saving a .csv and .npy information because the consequence.
Necessary! The VIIRS dataset contains a extremely beneficial variable cf_cvg, which describes the whole variety of observations that went into every pixel (cloud-free pixels). In essence, it’s a top quality flag. The larger this quantity, the upper high quality we’ve got. On this evaluation, when calculating the norm, we’ll filter out all of the pixels with cf_cvg≤1.
arrays, dates, rads = [], [], []
if information.measurement().getInfo()!=0:
data_list = information.toList(information.measurement())
for i in vary(data_list.measurement().getInfo()):
array, date = to_array(data_list,i, roi)rads.append(array['avg_rad'][np.where(array['cf_cvg']>1)].imply())
dates.append(date)
if date>=pd.to_datetime('2022-01-01'):
arrays.append(array['avg_rad'])
print(f'Index: {i+1}/{data_list.measurement().getInfo()+1}')
df = pd.DataFrame({'date': dates, 'avg_rad':rads})
np.save(f'{metropolis}.npy', arrays, allow_pickle=True)
df.to_csv(f'{metropolis}.csv', index=None)
The generated information of the format metropolis.csv with the avg_rad time collection inside are good for anomaly calculation. The method is very simple:
- Filter out observations earlier than 2022–02–01;
- Group by all of the observations by month (in complete — 12 teams);
- Take a imply;
- Subtract the imply from the observations after 2022–02–01 for every month respectively.
df = pd.read_csv(f'{metropolis}.csv')
df.date = pd.to_datetime(df.date)
ts_lon = df[df.date<pd.to_datetime('2022-01-01')].set_index('date')
means = ts_lon.groupby(ts_lon.index.month).imply()ts_short = df[df.date>=pd.to_datetime('2022-01-01')].set_index('date')
ts_short['month'] = ts_short.index.month
anomaly = ts_short['avg_rad']-ts_short['month'].map(means['avg_rad'])
Our final step to truly see the primary result’s constructing two subplots: a map + anomalies time collection. We’re not going to do any static maps at the moment. To implement a GIF, let’s construct a nested perform drawing our subplots:
def plot(metropolis, arrays, dates, rads):
def replace(body):
im1.set_data(arrays[frame])info_text = (
f"Date: {pd.to_datetime(dates[frame]).strftime(format='%Y-%m-%d')}n"
)
textual content.set_text(info_text)
ax[0].axis('off')
im2.set_data(dates[0:frame+1], rads[0:frame+1])
ax[1].relim()
return [im1, im2]
colours = [(0, 0, 0), (1, 1, 0)]
cmap_name = 'black_yellow'
black_yellow_cmap = LinearSegmentedColormap.from_list(cmap_name, colours)
llim = -1
fig, ax = plt.subplots(1,2,figsize=(12,8), frameon=False)
im1 = ax[0].imshow(arrays[0], vmax=10, cmap=black_yellow_cmap)
textual content = ax[0].textual content(20, 520, "", ha='left', fontsize=14, fontname='monospace', coloration='white')
im2, = ax[1].plot(dates[0], rads[0], marker='o',coloration='black', lw=2)
plt.xticks(rotation=45)
ax[1].axhline(0, lw=3, coloration='black')
ax[1].axhline(0, lw=1.5, ls='--', coloration='yellow')
ax[1].grid(False)
ax[1].spines[['right', 'top']].set_visible(False)
ax[1].set_xlabel('Date', fontsize=14, fontname='monospace')
ax[1].set_ylabel('Common DNB radiance', fontsize=14, fontname='monospace')
ax[1].set_ylim(llim, max(rads)+0.1)
ax[1].set_xlim(min(dates), max(dates))
ani = animation.FuncAnimation(fig, replace, frames=27, interval=40)
ani.save(f'{metropolis}.gif', fps=0.5, savefig_kwargs={'pad_inches':0, 'bbox_inches': 'tight'})
plt.present()
The code above is likely to be difficult to know at first look. Nevertheless it’s really fairly easy:
- Defining the replace perform. This perform is utilized by matplotlib perform FuncAnimation. The concept is that it passes (provides) new information to the present plot and returns a brand new determine (body). An inventory of frames is then reworked to a GIF file.
- Making a customized coloration map. This one is the best. I merely don’t like the colours of the built-in matplotlib cmaps for this venture. Since we’re working with gentle within the present evaluation, let’s use black and yellow.
- Constructing and formating the plots. Every little thing else is only a common map + line plot with labels, limits and ticks formatting. Nothing particular.
Let’s see what we’ve received right here: