Catapult Project - Australian Geoscience Datacube API

Perform band maths and produce a Normalised Difference Vegetation Index (NDVI) file.



In [11]:

    
from pprint import pprint
%matplotlib inline
from matplotlib import pyplot as plt
import xarray
import datacube.api



In [12]:

    
dc = datacube.api.API()



In [13]:

    
alos2 = dc.get_dataset(product='gamma0', platform='ALOS_2', 
                      y=(-42.55,-42.57), x=(147.55,147.57), 
                      variables=['hh_gamma0', 'hv_gamma0'])



In [14]:

    
s1a = dc.get_dataset(product='gamma0', platform='SENTINEL_1A', 
                      y=(-42.55,-42.57), x=(147.55,147.57), 
                      variables=['vh_gamma0', 'vv_gamma0'])



In [15]:

    
hh = alos2.hh_gamma0
hv = alos2.hv_gamma0
vh = s1a.vv_gamma0
vv = s1a.vv_gamma0

Select the first time index and plot the first timeslice (only one timeslice in this example).



In [16]:

    
hv.isel(time=0).plot.hist(bins=256, range=(0.0,1.0))









    Out[16]:





(array([  2.71500000e+03,   2.19400000e+03,   2.32700000e+03,
          2.57200000e+03,   2.58400000e+03,   2.50300000e+03,
          2.20000000e+03,   1.97800000e+03,   1.72200000e+03,
          1.35500000e+03,   1.17600000e+03,   1.00900000e+03,
          8.39000000e+02,   6.74000000e+02,   5.28000000e+02,
          4.55000000e+02,   3.96000000e+02,   3.11000000e+02,
          2.45000000e+02,   2.27000000e+02,   1.75000000e+02,
          1.55000000e+02,   1.23000000e+02,   9.80000000e+01,
          1.04000000e+02,   6.50000000e+01,   6.30000000e+01,
          5.60000000e+01,   5.20000000e+01,   3.20000000e+01,
          3.60000000e+01,   3.30000000e+01,   2.50000000e+01,
          1.40000000e+01,   2.00000000e+01,   1.40000000e+01,
          7.00000000e+00,   7.00000000e+00,   8.00000000e+00,
          7.00000000e+00,   1.30000000e+01,   1.00000000e+01,
          6.00000000e+00,   7.00000000e+00,   7.00000000e+00,
          2.00000000e+00,   2.00000000e+00,   3.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   2.00000000e+00,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   1.00000000e+00,
          1.00000000e+00,   1.00000000e+00,   0.00000000e+00,
          2.00000000e+00,   1.00000000e+00,   1.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          1.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   1.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00]),
 array([ 0.        ,  0.00390625,  0.0078125 ,  0.01171875,  0.015625  ,
         0.01953125,  0.0234375 ,  0.02734375,  0.03125   ,  0.03515625,
         0.0390625 ,  0.04296875,  0.046875  ,  0.05078125,  0.0546875 ,
         0.05859375,  0.0625    ,  0.06640625,  0.0703125 ,  0.07421875,
         0.078125  ,  0.08203125,  0.0859375 ,  0.08984375,  0.09375   ,
         0.09765625,  0.1015625 ,  0.10546875,  0.109375  ,  0.11328125,
         0.1171875 ,  0.12109375,  0.125     ,  0.12890625,  0.1328125 ,
         0.13671875,  0.140625  ,  0.14453125,  0.1484375 ,  0.15234375,
         0.15625   ,  0.16015625,  0.1640625 ,  0.16796875,  0.171875  ,
         0.17578125,  0.1796875 ,  0.18359375,  0.1875    ,  0.19140625,
         0.1953125 ,  0.19921875,  0.203125  ,  0.20703125,  0.2109375 ,
         0.21484375,  0.21875   ,  0.22265625,  0.2265625 ,  0.23046875,
         0.234375  ,  0.23828125,  0.2421875 ,  0.24609375,  0.25      ,
         0.25390625,  0.2578125 ,  0.26171875,  0.265625  ,  0.26953125,
         0.2734375 ,  0.27734375,  0.28125   ,  0.28515625,  0.2890625 ,
         0.29296875,  0.296875  ,  0.30078125,  0.3046875 ,  0.30859375,
         0.3125    ,  0.31640625,  0.3203125 ,  0.32421875,  0.328125  ,
         0.33203125,  0.3359375 ,  0.33984375,  0.34375   ,  0.34765625,
         0.3515625 ,  0.35546875,  0.359375  ,  0.36328125,  0.3671875 ,
         0.37109375,  0.375     ,  0.37890625,  0.3828125 ,  0.38671875,
         0.390625  ,  0.39453125,  0.3984375 ,  0.40234375,  0.40625   ,
         0.41015625,  0.4140625 ,  0.41796875,  0.421875  ,  0.42578125,
         0.4296875 ,  0.43359375,  0.4375    ,  0.44140625,  0.4453125 ,
         0.44921875,  0.453125  ,  0.45703125,  0.4609375 ,  0.46484375,
         0.46875   ,  0.47265625,  0.4765625 ,  0.48046875,  0.484375  ,
         0.48828125,  0.4921875 ,  0.49609375,  0.5       ,  0.50390625,
         0.5078125 ,  0.51171875,  0.515625  ,  0.51953125,  0.5234375 ,
         0.52734375,  0.53125   ,  0.53515625,  0.5390625 ,  0.54296875,
         0.546875  ,  0.55078125,  0.5546875 ,  0.55859375,  0.5625    ,
         0.56640625,  0.5703125 ,  0.57421875,  0.578125  ,  0.58203125,
         0.5859375 ,  0.58984375,  0.59375   ,  0.59765625,  0.6015625 ,
         0.60546875,  0.609375  ,  0.61328125,  0.6171875 ,  0.62109375,
         0.625     ,  0.62890625,  0.6328125 ,  0.63671875,  0.640625  ,
         0.64453125,  0.6484375 ,  0.65234375,  0.65625   ,  0.66015625,
         0.6640625 ,  0.66796875,  0.671875  ,  0.67578125,  0.6796875 ,
         0.68359375,  0.6875    ,  0.69140625,  0.6953125 ,  0.69921875,
         0.703125  ,  0.70703125,  0.7109375 ,  0.71484375,  0.71875   ,
         0.72265625,  0.7265625 ,  0.73046875,  0.734375  ,  0.73828125,
         0.7421875 ,  0.74609375,  0.75      ,  0.75390625,  0.7578125 ,
         0.76171875,  0.765625  ,  0.76953125,  0.7734375 ,  0.77734375,
         0.78125   ,  0.78515625,  0.7890625 ,  0.79296875,  0.796875  ,
         0.80078125,  0.8046875 ,  0.80859375,  0.8125    ,  0.81640625,
         0.8203125 ,  0.82421875,  0.828125  ,  0.83203125,  0.8359375 ,
         0.83984375,  0.84375   ,  0.84765625,  0.8515625 ,  0.85546875,
         0.859375  ,  0.86328125,  0.8671875 ,  0.87109375,  0.875     ,
         0.87890625,  0.8828125 ,  0.88671875,  0.890625  ,  0.89453125,
         0.8984375 ,  0.90234375,  0.90625   ,  0.91015625,  0.9140625 ,
         0.91796875,  0.921875  ,  0.92578125,  0.9296875 ,  0.93359375,
         0.9375    ,  0.94140625,  0.9453125 ,  0.94921875,  0.953125  ,
         0.95703125,  0.9609375 ,  0.96484375,  0.96875   ,  0.97265625,
         0.9765625 ,  0.98046875,  0.984375  ,  0.98828125,  0.9921875 ,
         0.99609375,  1.        ]),
 <a list of 256 Patch objects>)

Plot a band just for fun



In [49]:

    
vv.isel(time=0).plot.imshow(cmap="spectral", clim=(0.0, 0.5))









    Out[49]:





<matplotlib.image.AxesImage at 0x13624c41a90>

We can also select a range on the spatial dimensions.

Insert some rsgislib segmentation magic below

The xarray can be turned into a Dataset, which can then be saved across all of the timeseries to a NetCDF file.



In [44]:

    
ds_vv = vv[:100,:100,:100].to_dataset(name='vv')
ds_vv.to_netcdf('vv.nc')

Transposing dimensions

Merging two datasets



In [51]:

    
alos2_array = dc.get_dataset(product='gamma0', platform='ALOS_2', 
                      y=(-42.55,-42.57), x=(147.55,147.57), 
                      variables=['hh_gamma0', 'hv_gamma0'])
s1a_array = dc.get_dataset(product='gamma0', platform='SENTINEL_1A', 
                      y=(-42.55,-42.57), x=(147.55,147.57), 
                      variables=['vh_gamma0', 'vv_gamma0'])

sar_array = s1a_array.merge(alos2_array)
print (sar_array)









    



<xarray.Dataset>
Dimensions:    (time: 1, x: 156, y: 187)
Coordinates:
  * time       (time) datetime64[ns] 2016-03-01T23:59:59
  * y          (y) float64 -4.731e+06 -4.731e+06 -4.731e+06 -4.731e+06 ...
  * x          (x) float64 1.311e+06 1.311e+06 1.311e+06 1.311e+06 1.311e+06 ...
Data variables:
    vh_gamma0  (time, y, x) float32 0.0164571 0.0410914 0.0317888 0.0190076 ...
    vv_gamma0  (time, y, x) float32 0.0891735 0.0691831 0.128744 0.122255 ...
    crs        int32 0
    hh_gamma0  (time, y, x) float32 0.0655217 0.143631 0.11165 0.116711 ...
    hv_gamma0  (time, y, x) float32 0.0378705 0.0161559 0.0106346 0.00677407 ...
Attributes:
    license: Creative Commons Attribution 4.0 International CC BY 4.0
    source: This data is a reprojection and retile of Landsat surface reflectance data from the USGS
    summary: These files are experimental, short lived, and the format will change.
    title: Experimental Data files From the Australian Geoscience Data Cube - DO NOT USE
    product_version: 0.0.0

We will clip and scale the image to improve the contrast of the visible bands.



In [61]:

    
fake_saturation = 0.05
clipped_sar = sar_array.where(sar_array<fake_saturation).fillna(fake_saturation)
max_val = clipped_sar.max(['y', 'x'])
scaled = (clipped_sar / max_val)



In [63]:

    
rgb = scaled.transpose('time', 'y', 'x')
rgb.dims









    Out[63]:





Frozen(SortedKeysDict({'y': 187, 'x': 156, 'time': 1}))



In [64]:

    
plt.imshow(rgb.isel(time=0))









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-64-292cc5735cee> in <module>()
----> 1 plt.imshow(rgb.isel(time=0))

C:\Users\Simon\Anaconda3\envs\py35\lib\site-packages\matplotlib\pyplot.py in imshow(X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, shape, filternorm, filterrad, imlim, resample, url, hold, data, **kwargs)
   3020                         filternorm=filternorm, filterrad=filterrad,
   3021                         imlim=imlim, resample=resample, url=url, data=data,
-> 3022                         **kwargs)
   3023     finally:
   3024         ax.hold(washold)

C:\Users\Simon\Anaconda3\envs\py35\lib\site-packages\matplotlib\__init__.py in inner(ax, *args, **kwargs)
   1809                     warnings.warn(msg % (label_namer, func.__name__),
   1810                                   RuntimeWarning, stacklevel=2)
-> 1811             return func(ax, *args, **kwargs)
   1812         pre_doc = inner.__doc__
   1813         if pre_doc is None:

C:\Users\Simon\Anaconda3\envs\py35\lib\site-packages\matplotlib\axes\_axes.py in imshow(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, shape, filternorm, filterrad, imlim, resample, url, **kwargs)
   4945                               resample=resample, **kwargs)
   4946 
-> 4947         im.set_data(X)
   4948         im.set_alpha(alpha)
   4949         if im.get_clip_path() is None:

C:\Users\Simon\Anaconda3\envs\py35\lib\site-packages\matplotlib\image.py in set_data(self, A)
    447         if (self._A.dtype != np.uint8 and
    448                 not np.can_cast(self._A.dtype, np.float)):
--> 449             raise TypeError("Image data can not convert to float")
    450 
    451         if (self._A.ndim not in (2, 3) or

TypeError: Image data can not convert to float



In [ ]:

    
import matplotlib.image
matplotlib.image.imsave('vv_hh_vh.png', rgb.isel(time=16))

Behind the scenes

The ndvi result is performed by chaining a series of operations together on a per-chunk basis. The execution tree, including the masking on the nodata value done by the API, ends up being the same for each chunk. The graph can be read from the bottom up, with the ouput array chunks at the top. (Double-click the tree-graph image below to zoom)



In [ ]:

    
sar.data.visualize()

If we look at a single chunk, the NDVI calculation can be seen where the lines cross over to the add and sub circles. Some optimizarion has taken place: the div operation has been combined with another inline function, and the other chunks have been discarded.



In [ ]:

    
partial = vv[0,0,0]
partial.data.visualize(optimize_graph=True)