In this notebook you will learn more about the Kyrgyzstan Data Cube (KDC) products and how to use geopandas to analyse or view the region that you want.
This notebook will require a number of external Python modules to run, which will be imported at the top of the notebook as his convention.
import datetime
import sys
import geopandas
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from datacube import Datacube
from rasterio.features import geometry_mask
sys.path.append('/eOsphere/bin')
import dcFunctions
As we have been showed in other notebooks, we need to initialise a Data Cube in order to have access to the KDC, and this can be done as follows:
dc = Datacube()
We use geopandas to load the Raions shapefiles; inside there will be a poligon for each raion.
Use the command print(sf) to see what's inside the shapefile.
sf = geopandas.read_file('data/kgz_admbnda_adm2_moes_20181119.shp', encoding='utf-8')
sf
As showed in previous notebooks, you can select a specific Raion using the index of the geopandas dataframe.
# Select a Raion
region_id = 15
# Get the coordinate of the box of the Raion selected
minX, minY, maxX, maxY = sf.bounds.loc[region_id]
Define query to load data from the Data Cube.
In this case, we are loading the already made 10day median NDVI, in a raion (region_id=15), for the last 10day of June.
start_date, end_date = '2019-06-21','2019-06-30'
query = {'x': [minX,maxX],
'y': [minY,maxY],
'time': [start_date, end_date],
'measurements': ['ndvi'],
'output_crs': 'EPSG:32643',
'resolution': [10,10],
'dask_chunks': {'time': 1, 'x': 1024, 'y': 1024}
}
ds = dc.load(product='indices', **query)
crs = ds.crs
shape_ds = ds.geobox.shape
affine_ds = ds.geobox.affine
Let's scale the data so to have the NDVI from -1 to +1.
print(ds)
ds = ds.isel(time=0)
ds = dcFunctions.scal_2_tru(ds)
Here we change the CRS of the shapefile to make sure that the shapefile and the data are in the same projection.
sf_utm = sf.to_crs(epsg=32643)
Now we plot the NDVI with the borders of the raion selected.
fig, ax = plt.subplots()
ds.ndvi.plot(ax=ax, cmap='RdYlGn')
sf_utm.plot(ax=ax, facecolor='none', edgecolor='black')
# This is just to make sure that we are using the same projection of the data
gdf = sf.iloc[[region_id]].to_crs(crs)
# Create a mask from the shapefile
mask = geometry_mask(
gdf.geometry,
out_shape = shape_ds,
transform = affine_ds,
all_touched = False,
invert = False)
mask = xr.DataArray(mask, dims=('y', 'x'))
ds = ds.where(~mask)
Plot the cutted data with the shapefile of the raions
fig, ax = plt.subplots()
ds.ndvi.plot(ax=ax, cmap='RdYlGn')
sf_utm.plot(ax=ax, facecolor='none', edgecolor='black')
Plot the cutted data with the shapefile of the selected raion
fig, ax = plt.subplots()
ds.ndvi.plot(ax=ax, cmap='RdYlGn')
gdf.plot(ax=ax, facecolor='none', edgecolor='black')
In the following cells we will use MODIS data to show how the data can be cutted using multiple polygons.
# Load the shapefile and find the bounds
sf = geopandas.read_file('data/kgz_admbnda_adm2_moes_20181119.shp')
minX, minY = sf.bounds['minx'].min(), sf.bounds['miny'].min()
maxX, maxY = sf.bounds['maxx'].max(), sf.bounds['maxy'].max()
# Query to load MODIS data
start_date, end_date = '2019-06-21','2019-06-30'
query = {'x': [minX,maxX],
'y': [minY,maxY],
'time': [start_date, end_date],
'measurements': ['red','green','blue'],
'output_crs': 'EPSG:32643',
'resolution': [-500,500],
'dask_chunks': {'time': 1, 'x': 1024, 'y': 1024}
}
ds = dc.load(product='modis_nbar', **query)
# Save some parameters for later on
crs = ds.crs
shape_ds = ds.geobox.shape
affine_ds = ds.geobox.affine
print(ds)
# Remove nodata value
ds = ds.where(ds!=32767)
# Compute the median for the 3 channels/bands
ds = ds.median(dim='time').compute()
# Apply a correction (for the colours)
ds = ds**0.5
print(ds)
RGB of the whole country with the shapefile of the Raions
fig, ax = plt.subplots(figsize = (15,7))
ds[['red','green','blue']].to_array().plot.imshow(ax=ax, robust=True)
sf_utm.plot(ax=ax, facecolor='none', edgecolor='red')
Cut the image selecting only two Raions (15 and 33 are the index in the geopandas dataframe)
# List all the raions that you want
region_ids = [15, 33]
mask = np.ones(shape_ds, dtype=bool)
for region in region_ids:
# This is just to make sure that we are using the same projection of the data
gdf = sf.iloc[[region]].to_crs(crs)
# Create a mask from the shape file
temp = geometry_mask(
gdf.geometry,
out_shape = shape_ds,
transform = affine_ds,
all_touched = False,
invert = False)
mask = mask & temp
mask = xr.DataArray(mask, dims=('y', 'x'))
# Cropping our image
ds = ds.where(~mask)
Plotting the cutted image
fig, ax = plt.subplots(figsize = (15,7))
ds[['red','green','blue']].to_array().plot.imshow(ax=ax, robust=True)
sf_utm.plot(ax=ax, facecolor='none', edgecolor='red')
The same can be applied for the Sentinel data. We can load from the RGB product
dc.list_measurements().loc['rgb']
Load the raion shapefile and selecting a specific soum
# Select a raion
sf = geopandas.read_file('data/kgz_admbnda_adm2_moes_20181119.shp')
region_id = 15
minX, minY, maxX, maxY = sf.bounds.loc[region_id]
Setting the query for the Data Cube to load data from high resolution data
start_date, end_date = '2019-06-01','2019-06-30'
query = {'x': [minX,maxX],
'y': [minY,maxY],
'time': [start_date, end_date],
'measurements': ['red','green','blue'],
'output_crs': 'EPSG:32643',
'resolution': [10,10],
'dask_chunks': {'time': 1, 'x': 1024, 'y': 1024}
}
Load the data and make sure that our raion (gdf) has the same projection of our data (crs = EPSG:32643)
ds = dc.load(product='rgb', **query)
print(ds)
crs = ds.crs
gdf = sf.iloc[[region_id]].to_crs(crs)
Select the only one image/time snapshot available in this case.
ds = ds.isel(time=0)
ds = dcFunctions.scal_2_tru(ds)
print(ds)
Plotting the data and the boundary of the selected raion
fig, ax = plt.subplots(figsize = (15,7))
ds[['red','green','blue']].to_array().plot.imshow(ax=ax, robust=True)
gdf.plot(ax=ax, facecolor='none', edgecolor='red')
This notebook will attempt to show some of the long term time analysis processes that can be run using the Data Cube. The Data Cube is designed to allow for easy analysis of large volumes of Analysis Ready Data in both spatial and temporal dimensions. In this notebook, we will look at using Sentinel-2 data to analyse the Modified Soil Adjusted Vegetation Index (MSAVI) for a group of fields near Bishkek over time as well as find how the Normalised Difference Vegetation Index (NDVI) and Normalised Difference Water Index (NDWI) changes over the entire Naryn oblast in Kyrgyzstan for a whole year.
To start off with, we will need to load in the modules required to do this analysis, which can be done below:
import sys
import datacube
import numpy as np
import matplotlib.pyplot as plt
sys.path.append('/eOsphere/bin')
import dcFunctions
The next thing to do is consider the area which will be analysed. We will be looking at some fields located to the West of Bishkek. Its latitude range is between 42.90 and 42.91 and its longitude range is between 74.39 and 74.91. Below is a picture of the area containing the fields, which is constrained by the red box.
For the purposes of this exercise, the bounding box we will use is the latitude and longitude extent range.
Now we have enough information to initialise a Data Cube instance, and create a query for the search. Try to do this below, with an answer for how to do this in the hidden solutions box.
dc = datacube.Datacube()
query = {'x': [74.386, 74.391],
'y': [42.9, 42.91],
'measurements': ['red', 'nir', 'mask'],
'output_crs': 'EPSG:32643',
'resolution': [10, 10]
}
Once the query has been done, it's time to load in the data. In this analysis, we will be using Sentinel-2 10m data. As we will be looking at the data for the entire Sentinel-2 archive currently in the Kyrgyz Data Cube, it is not necessary to set a time range, and so the data can be loaded in as shown below:
ds = dc.load(product='s2_10m', **query)
ds = ds.sortby('time')
print(ds)
This data must be masked out, so that anything that isn't land or snow is ignored. This can be done the same way as shown in previous notebooks, bearing in mind that the values for land and snow in the mask DataArray are 5 and 7 respectively. Try masking the dataset below, looking at the solution in the box if you need help.
ds = ds.where((ds.mask == 5) | (ds.mask == 7))
To find the median average pasture for this farmland area, we will need to take the median across the 'x' and 'y' dimensions which can be done in the following way:
ds = ds.median(['x', 'y'])
print(ds)
Now that the exact coordinate of Station 246 has been identified in the data and the data has been cloud masked, it is time to create the MSAVI2 index. This index is similar to NDVI, in that it takes into account both the red and near-infrared bands from satellite imagery, but it accounts better for areas with large amounts of bare soil. The equation for implementing the MSAVI was derived by Qi et al. (1994) and is as follows:
$MSAVI = 2 \ \rho_{NIR} + 1 - \frac{\sqrt{(2 \ \rho_{NIR} + 1)^2 - 8 \ (\rho_{NIR} - \rho_{RED})}}{2} $
where $\rho_{NIR}$ is the near-infrared band surface reflectance and $\rho_{RED}$ is the red band surface reflectance. This equation can be calculated using Python as shown below:
msavi = ((2 * ds.nir) + 1 - ((2 * ds.nir + 1)**2 - 8 * (ds.nir - ds.red))**0.5) / 2
This has created an xarray DataArray of MSAVI points. These points can be plotted using Matplotlib and the .scatter() function, using the the points created in the MSAVI array and lining them up with the time array given by ds.time.
plt.scatter(ds.time, msavi, marker='.')
Now we can do some timeseries averaging to find useful information for a larger area over the course of a year. In this instance, we will plot the average value for all of Naryn oblast for a number of different indices over the course of a year. This will be done using composite Sentinel-2/Landsat8 data at 100m resolution. This will require a change to the query used in the previous part of the notebook.
Naryn has the following bbox coordinates:
It is also important to remember that the pixel resolution is now 100m. This results in the following query:
query = {'x': [73.7128139, 77.9073308],
'y': [40.2809711, 42.4530034],
'time': ['2019-01-01', '2019-12-31'],
'measurements': ['ndvi', 'ndwi'],
'output_crs': 'EPSG:32643',
'resolution': [100, 100]
}
This query can then be used to get all the 10 Day composite data in for Zavkhan in 2019. The product required to do this is modis_indices_250m and otherwise it is done in the usual way. Try to do this below, looking at the hidden solution if required.
ds = dc.load(product='month_indices', **query)
print(ds)
The loaded in data must now be scaled from the uint16 values it is stored as to the familiar float values that are expected when for indices such as NDVI, NDWI and NDDI. This can be done using the function scal_2_tru in dcFunctions, as shown below:
ds = dcFunctions.scal_2_tru(ds)
Next we must find the median value for NDVI and NDWI at each time step. This can be done by once again using the .median() method on the dataset itself, but instead of specifying the dimension that the median will be taken over as the "time" dimension, we will instead specify this using just square brackets, to show that we want the median over both the "x" and "y" dimensions. This is done as shown below:
ds = ds.median(['y','x'])
print(ds)
Next we can plot the average NDVI and NDWI over the course of a year for Naryn oblast. This is done using matplotlib, and to plot both lines on the same graph, a little more setup must be done compared to previous examples.
To do this we have to set up the plot figure before initialising the plot, using the plt.subplots() function. For each plot we specify the axes to plot the figure on, to ensure both plots overlap on the same axes. We can also use the "label" function, which when combined with plt.legend() creates an easy label box for both lines.
The below example illustrates just some of the things that can be customised in plotting with matplotlib, although it is possible to do a very wide variety of things indeed.
fig, ax = plt.subplots(1, 1, figsize=[8,5])
ds['ndvi'].plot.line(ax=ax, label='NDVI')
ds['ndwi'].plot.line(ax=ax, label='NDWI')
plt.legend(loc='upper left')
plt.show()