In [1]:
import datacube.api
import numpy
import fiona
import geopandas
import rasterio
from affine import Affine
from rasterio.crs import from_string
import osr
from skimage import exposure
from image_processing.segmentation import Segments, rasterise_vector
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
dc = datacube.api.API()
Open the vector file and get the bounding box
In [3]:
vfname = '/g/data2/v10/public/water-example/sample-water-bodies.shp'
src = fiona.open(vfname, 'r')
xidx = (src.bounds[0], src.bounds[2])
yidx = (src.bounds[-1], src.bounds[1])
gdf = geopandas.read_file(vfname)
gdf.plot()
Out[3]:
Using the co-ordinates of the vector file's bunding box, query the datacube.
In [4]:
nbar = dc.get_dataset(product='nbar', platform='LANDSAT_8',
y=yidx, x=xidx)
nbar
Out[4]:
Extract and convert the co-ordinate reference system into rasterio's preferred format
In [5]:
sr = osr.SpatialReference()
sr.ImportFromWkt(nbar.crs.spatial_ref)
crs = from_string(sr.ExportToProj4())
print crs
Create an affine/geotransformtion matrix.
We make an assumption that the spatial array is uniformally spaced in order for this to work.
In [6]:
pix_x = nbar.x.values[1] - nbar.x.values[0]
pix_y = nbar.y.values[1] - nbar.y.values[0]
ulx = nbar.x.values[0] - pix_x / 2.0
uly = nbar.y.values[0] - pix_y /2.0
transform = Affine.from_gdal(*[ulx, pix_x, 0, uly, 0, pix_y])
print transform
Determine the pixel size in hectares. We'll use this later on.
In [7]:
ha = pix_x **2 / 10000.0
print "Pixel area in hectares: {}".format(ha)
In [8]:
dims = nbar.band_4.shape
In [9]:
img = numpy.zeros((3, dims[1], dims[2]), dtype='float32')
img[0] = nbar.band_5[28]
img[1] = nbar.band_4[28]
img[2] = nbar.band_3[28]
img[img == -999] = numpy.nan
img /= 10000
In [10]:
scl = exposure.equalize_hist(img, mask=numpy.isfinite(img))
In [11]:
plt.imshow(scl.transpose(1, 2, 0))
Out[11]:
In [12]:
ras = rasterise_vector(vfname, shape=dims[1:], crs=crs, transform=transform)
seg = Segments(ras)
print "Number of segments: {}".format(seg.n_segments)
For the NIR, Green, & Red bands, calculate statistics for every segment.
In [13]:
dat = nbar.band_5[28].values.astype('float32')
dat[dat == -999] = numpy.nan
nir_stats = seg.basic_statistics(dat, nan=True, scale_factor=ha)
dat = nbar.band_4[28].values.astype('float32')
dat[dat == -999] = numpy.nan
red_stats = seg.basic_statistics(dat, nan=True, scale_factor=ha)
dat = nbar.band_3[28].values.astype('float32')
dat[dat == -999] = numpy.nan
green_stats = seg.basic_statistics(dat, nan=True, scale_factor=ha)
In [14]:
nir_stats.head(10)
Out[14]:
In [15]:
bboxes = seg.bounding_box()
In [16]:
sid = 185
window = bboxes[sid]
ys, ye = window[0]
xs, xe = window[1]
print window
In [17]:
subs = img[:, ys:ye, xs:xe]
scl_subs = exposure.equalize_hist(subs, mask=numpy.isfinite(subs))
plt.title(src[sid - 1]['properties']['FEATURETYP'] + ' (???)')
plt.imshow(scl_subs.transpose(1, 2, 0))
Out[17]:
Find a segment with a low mean value.
In [18]:
nir_stats[nir_stats['Mean'] < 500]
Out[18]:
In [19]:
sid = 357
window = bboxes[sid]
ys, ye = window[0]
xs, xe = window[1]
print window
In [20]:
subs = img[:, ys:ye, xs:xe]
scl_subs = exposure.equalize_hist(subs, mask=numpy.isfinite(subs))
plt.title(src[sid - 1]['properties']['FEATURETYP'])
plt.imshow(scl_subs.transpose(1, 2, 0))
Out[20]:
In [21]:
sid = 262
window = bboxes[sid]
ys, ye = window[0]
xs, xe = window[1]
print window
In [22]:
subs = img[:, ys:ye, xs:xe]
scl_subs = exposure.equalize_hist(subs, mask=numpy.isfinite(subs))
plt.title(src[sid - 1]['properties']['NAME'])
plt.imshow(scl_subs.transpose(1, 2, 0))
Out[22]:
In [23]:
dat = nbar.band_6[0].values.astype('float32')
dat[dat == -999] = numpy.nan
swir_stats = seg.basic_statistics(dat, nan=True, scale_factor=ha)
swir_stats['timestamp'] = nbar.time[0].values
for i in range(1, dims[0]):
dat = nbar.band_6[i].values.astype('float32')
dat[dat == -999] = numpy.nan
stats = seg.basic_statistics(dat, nan=True, scale_factor=ha)
stats['timestamp'] = nbar.time[i].values
swir_stats = swir_stats.append(stats)
In [24]:
swir_stats.set_index('timestamp', inplace=True)
In [25]:
print "Number of records: {}".format(swir_stats.shape[0])
swir_stats.head(10)
Out[25]:
In [26]:
sid = 262
roi = swir_stats[swir_stats['Segment_IDs'] == sid]
roi = roi[numpy.isfinite(roi['Mean'])]
roi['Mean'].plot(title=src[sid - 1]['properties']['NAME'])
Out[26]:
In [27]:
sid = 357
roi = swir_stats[swir_stats['Segment_IDs'] == sid]
roi = roi[numpy.isfinite(roi['Mean'])]
roi['Mean'].plot(title=src[sid - 1]['properties']['FEATURETYP'])
Out[27]:
In [28]:
sid = 185
roi = swir_stats[swir_stats['Segment_IDs'] == sid]
roi = roi[numpy.isfinite(roi['Mean'])]
roi['Mean'].plot(title=src[sid - 1]['properties']['FEATURETYP'])
Out[28]:
In [29]:
sid = 262
window = bboxes[sid]
ys, ye = window[0]
xs, xe = window[1]
roi = swir_stats[swir_stats['Segment_IDs'] == sid]
roi = roi[numpy.isfinite(roi['Mean'])]
In [30]:
wh1 = (roi['Mean'] < 1100) & (roi['Total'] > 110000000)
roi[wh1]
Out[30]:
In [31]:
wh2 = (roi['Mean'] > 2900) & (roi['StdDev'] < 400)
roi[wh2]
Out[31]:
In [32]:
fig, axes = plt.subplots(ncols=2)
timestamp = roi[wh1].index[0]
subs[0] = nbar.band_5.loc[timestamp][ys:ye, xs:xe]
subs[1] = nbar.band_4.loc[timestamp][ys:ye, xs:xe]
subs[1] = nbar.band_3.loc[timestamp][ys:ye, xs:xe]
scl_subs = exposure.equalize_hist(subs, mask=subs != -999)
axes[0].set_title(timestamp)
axes[0].imshow(scl_subs.transpose(1, 2, 0))
timestamp = roi[wh2].index[0]
subs[0] = nbar.band_5.loc[timestamp][ys:ye, xs:xe]
subs[1] = nbar.band_4.loc[timestamp][ys:ye, xs:xe]
subs[1] = nbar.band_3.loc[timestamp][ys:ye, xs:xe]
scl_subs = exposure.equalize_hist(subs, mask=subs != -999)
axes[1].set_title(timestamp)
axes[1].imshow(scl_subs.transpose(1, 2, 0))
Out[32]:
In [ ]: