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

OBJECTIDint64

1 - 37

shapeNameobject

Andaman & Nicobar Islands2.7%

Andhra Pradesh2.7%

35 others94.6%

Andaman & Nicobar Islands

Andhra Pradesh

Arunachal Pradesh

Assam

Bihar

Chandigarh

Chhattisgarh

Dadra & Nagar Haveli

Daman & Diu

Delhi

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

/opt/conda/lib/python3.9/site-packages/pyproj/crs/crs.py:131: FutureWarning:

'+init=<authority>:<code>' syntax is deprecated. '<authority>:<code>' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6

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