Introduction to remote sensing: From satellite images to tabular data
Setup
import numpy
import matplotlib.pyplot as plt
import geopandas as gpd
import contextily as cx
import plotly.express as px
import plotly.graph_objects as go
import rasterio
from rasterio import plot as rioplot
from rasterio.mask import mask
import rasterstats
Load vector data
polygonsFile = gpd.read_file('/work/vector-data/india37map.geojson')
polygonsFile.plot();
polygonsFile.columns
polygonsFile.crs
polygonsFile_bbox = polygonsFile.total_bounds
Load raster data
rasterFile = rasterio.open('/work/raster-data/NTL_2020.tif')
rasterFile_window = rasterFile.window(*polygonsFile_bbox)
rasterFile_clipped = rasterFile.read(1, window = rasterFile_window)
rasterFile.crs
rasterFile.shape
rasterFile.bounds
rasterFile_extent = numpy.asarray(rasterFile.bounds)[[0,2,1,3]]
plt.figure(figsize=(10,10))
plt.imshow(rasterFile_clipped, extent=polygonsFile_bbox[[0,2,1,3]], cmap='magma')
polygonsFile.boundary.plot(ax=plt.gca(), color='skyblue', linewidth = 0.2);
# Recommended cmaps: hot, inferno, plasma, magma, cividis
Compute zonal statistics
polygonsFile_transform = rasterFile.window_transform(rasterFile_window)
def clean_mask(geom, dataset=rasterFile, **mask_kw):
mask_kw.setdefault('crop', True)
mask_kw.setdefault('all_touched', True)
mask_kw.setdefault('filled', False)
masked, mask_transform = mask(dataset=dataset, shapes=(geom,),
**mask_kw)
return masked
polygonsFile['rasterValues_sum'] = polygonsFile.geometry.apply(clean_mask).apply(numpy.ma.sum)
polygonsFile['rasterValues_median'] = polygonsFile.geometry.apply(clean_mask).apply(numpy.ma.median)
polygonsFile
Plot zonal statistics
gdf = gpd.read_file("/work/vector-data/india37map-OPT.geojson")
gdf['rasterValues_sum'] = polygonsFile['rasterValues_sum']
gdf['rasterValues_median'] = polygonsFile['rasterValues_median']
gdf.columns
px.histogram(gdf, x = 'rasterValues_sum', nbins = 16, hover_name= 'shapeName', marginal='rug')
gdf.explore(
column='rasterValues_sum',
tooltip=['shapeName', 'rasterValues_sum'],
scheme = 'FisherJenks', # Quantiles, EqualInterval, BoxPlot, FisherJenks
cmap = 'inferno', # hot, cividis, plasma, magma, inferno
legend=True,
tiles = 'CartoDB dark_matter', # OpenStreetMap, Stamen Terrain, Stamen Toner, Stamen Watercolor, CartoDB positron, CartoDB dark_matter
style_kwds =dict(color="gray", weight=0.2),
legend_kwds=dict(colorbar=False)
)
px.histogram(gdf, x = 'rasterValues_median', nbins = 16, hover_name= 'shapeName', marginal='rug')
gdf.explore(
column='rasterValues_median',
tooltip=['shapeName', 'rasterValues_median'],
scheme = 'FisherJenks', # Quantiles, EqualInterval, BoxPlot, FisherJenks
cmap = 'inferno', # hot, cividis, plasma, magma, inferno
legend=True,
tiles = 'CartoDB dark_matter', # OpenStreetMap, Stamen Terrain, Stamen Toner, Stamen Watercolor, CartoDB positron, CartoDB dark_matter
style_kwds =dict(color="gray", weight=0.2),
legend_kwds=dict(colorbar=False)
)
Plot SOL per area
polygonsFile = polygonsFile.to_crs({'init': 'epsg:24378'})
polygonsFile.crs
polygonsFile['areaKM2'] = polygonsFile['geometry'].area/ 10**6
polygonsFile['SOLpa'] = polygonsFile['rasterValues_sum'] / polygonsFile['areaKM2']
gdf['SOLpa'] = polygonsFile['SOLpa']
px.histogram(gdf, x = 'SOLpa', nbins = 16, hover_name= 'shapeName', marginal='rug')
gdf.explore(
column='SOLpa',
tooltip=['shapeName', 'rasterValues_median'],
scheme = 'FisherJenks', # Quantiles, EqualInterval, BoxPlot, FisherJenks
cmap = 'inferno', # hot, cividis, plasma, magma, inferno
legend=True,
tiles = 'CartoDB dark_matter', # OpenStreetMap, Stamen Terrain, Stamen Toner, Stamen Watercolor, CartoDB positron, CartoDB dark_matter
style_kwds =dict(color="gray", weight=0.2),
legend_kwds=dict(colorbar=False)
)
Plot log (1+SOLpa)
polygonsFile['ln_SOLpa'] = numpy.log(1 + polygonsFile['SOLpa'])
gdf['ln_SOLpa'] = polygonsFile['ln_SOLpa']
px.histogram(gdf, x = 'ln_SOLpa', nbins = 16, hover_name= 'shapeName', marginal='rug')
gdf.explore(
column='ln_SOLpa',
tooltip=['shapeName', 'rasterValues_median'],
scheme = 'FisherJenks', # Quantiles, EqualInterval, BoxPlot, FisherJenks
cmap = 'inferno', # hot, cividis, plasma, magma, inferno
legend=True,
tiles = 'CartoDB dark_matter', # OpenStreetMap, Stamen Terrain, Stamen Toner, Stamen Watercolor, CartoDB positron, CartoDB dark_matter
style_kwds =dict(color="gray", weight=0.2),
legend_kwds=dict(colorbar=False)
)
Export tabular data
polygonsFile
polygonsFile.drop('geometry',axis=1).to_csv('tabular-data/NTL.csv')