How much green is Delhi as on 13 Oct 2022?

Introduction

The India State of Forests Report (ISFR) 2021 [1], showed that the green cover of Delhi has increased from 21.88% in 2019 to 23.06% in 2021, of its geographical area as observed from satellite imageries of Delhi [2]. This was a welcome news for the city struggling with severe pollution and rising population, which makes it necessary to monitor the city's green cover and keep the city liveable.

This sample shows the capabilities of spectral indices such as Normalized Difference Vegetation index (NDVI) for the calculation of green cover in Delhi on 21 October 2022 using Landsat 8 imagery.

Necessary Imports

%matplotlib inline

import pandas as pd
from datetime import datetime
from IPython.display import Image
from IPython.display import HTML
import matplotlib.pyplot as plt

import arcgis
from arcgis.gis import GIS
from arcgis.raster.functions import apply, clip, remap, colormap
from arcgis.geocoding import geocode

Connect to your GIS

gis = GIS("home")

Get the data for analysis

Here we're getting the multispectral landsat imagery item in ArcGIS Online.

landsat_item = gis.content.get('d9b466d6a9e647ce8d1dd5fe12eb434b')
landsat = landsat_item.layers[0]
landsat_item
Multispectral Landsat
Landsat multispectral and multitemporal imagery with on-the-fly renderings and indices for visualization and analysis. The Landsat 8 and 9 imagery in this layer is updated daily and is directly sourced from the USGS Landsat collection on AWS.Imagery Layer by esri
Last Modified: July 01, 2022
3 comments, 855238 views

Search for India State Boundaries 2022 layer in ArcGIS Online. This layer has all the state boundaries for India. The boundary of Delhi can be filtered from the layer, as this notebook focuses on the city's green cover.

boundaries = gis.content.get('b66401aa25074b098aaa571b10c1b21f')
state_boundaries = boundaries.layers[1]
area = geocode("New Delhi, India", out_sr=landsat.properties.spatialReference)[0]
landsat.extent = area['extent']

Extracting Landsat imagery for Delhi Region

In State Boundary layer, OBJECTID for Delhi is 7 which is used below. Also, it is important to add extent to the geometry of selected boundary.

delhi = state_boundaries.query(where='OBJECTID=7')
delhi_geom = delhi.features[0].geometry
delhi_geom['spatialReference'] = {'wkid':3857}
delhi.features[0].extent = area['extent']

Filter imageries based on cloud cover and Acquisition Date

In order to have good result, it is important to select cloud free imagery from the image collection for a specified time duration. In this example we have selected all the imageries captured between 1 October, 2022 to 31 December, 2022 with cloud cover less than or equal to 5% for Delhi.

from datetime import datetime

selected = landsat.filter_by(where="(Category = 1) AND (cloudcover <=0.05)",
                             time=[datetime(2022, 10, 1), datetime(2022, 12, 31)],
                             geometry=arcgis.geometry.filters.intersects(delhi_geom))

df = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover", 
                    order_by_fields="AcquisitionDate").sdf
df['AcquisitionDate'] = pd.to_datetime(df['AcquisitionDate'], unit='ms')
df
OBJECTIDAcquisitionDateGroupNameCloudCoverSHAPE
038501172022-10-04 05:25:39LC08_L1TP_147040_20221004_20221012_02_T1_MTL0.0031{"rings": [[[8624803.6956, 3442830.7743000016]...
138624012022-10-13 05:19:21LC08_L1TP_146040_20221013_20221020_02_T1_MTL0.0013{"rings": [[[8793415.397599999, 3411368.579099...
238809762022-10-20 05:25:36LC08_L1TP_147040_20221020_20221101_02_T1_MTL0.0007{"rings": [[[8626312.848000001, 3442797.9736],...
338667512022-10-21 05:19:18LC09_L1TP_146040_20221021_20221021_02_T1_MTL0.0137{"rings": [[[8775485.517, 3347563.3280000016],...
438739562022-10-28 05:25:29LC09_L1TP_147040_20221028_20221028_02_T1_MTL0.023{"rings": [[[8626123.000799999, 3442751.5911],...
538916542022-10-29 05:19:29LC08_L1TP_146040_20221029_20221108_02_T1_MTL0.0247{"rings": [[[8796077.917599998, 3442767.058200...
639369482022-11-22 05:19:19LC09_L1TP_146040_20221122_20221122_02_T1_MTL0.0024{"rings": [[[8786887.706999999, 3392825.541599...
739574712022-11-30 05:19:25LC08_L1TP_146040_20221130_20221206_02_T1_MTL0.0004{"rings": [[[8796814.0528, 3442741.770499997],...
839677102022-12-07 05:25:35LC08_L1TP_147040_20221207_20221213_02_T1_MTL0.0017{"rings": [[[8624721.7905, 3442832.8721999973]...
939607782022-12-08 05:19:18LC09_L1TP_146040_20221208_20221208_02_T1_MTL0.0097{"rings": [[[8798761.3933, 3442669.745099999],...
1039705252022-12-15 05:25:29LC09_L1TP_147040_20221215_20221215_02_T1_MTL0.0009{"rings": [[[8626467.222199999, 3442743.814499...
1139858402022-12-16 05:19:19LC08_L1TP_146040_20221216_20221227_02_T1_MTL0.0271{"rings": [[[8797672.5849, 3442843.060800001],...
1239820372022-12-24 05:19:15LC09_L1TP_146040_20221224_20221224_02_T1_MTL0.0135{"rings": [[[8799790.6417, 3442584.7436999977]...

Selecting imagery dated 21 October, 2022 since the OBJECTID value in the query below is for that date.

delhi_image = landsat.filter_by('OBJECTID=3866751') # 2022-10-21

Applying Natural color to verify the quality of image.

apply(delhi_image, 'Natural Color with DRA')
<ImageryLayer url:"https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer">

Creating NDVI composite for the filtered imagery

In the Landsat layer properties, pre defined "NDVI Raw" function is applied to get the NDVI composite.

ndvi_colorized = apply(delhi_image, 'NDVI Raw')
ndvi_colorized
<ImageryLayer url:"https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer">

Clipping the NDVI composite for Delhi and setting the extent of the generated raster.

delhi_clip = clip(raster=ndvi_colorized, geometry= delhi_geom)
delhi_clip.extent = area['extent']
delhi_clip
<ImageryLayer url:"https://landsat2.arcgis.com/arcgis/rest/services/Landsat/MS/ImageServer">

Masking NDVI composite for green area calculation

"remap" function is used to define the NDVI range for agricultural land and forest. The NDVI values between 0.4 - 0.5 represents agricultural land whereas NDVI values between 0.5 - 1 shows forest/ tree cover [2].

threshold_val = 0.5
masked = colormap(remap(delhi_clip, 
                        input_ranges=[0.4,threshold_val,     # agricultural land
                                     threshold_val, 1],      # forest area/ tree cover
                        output_values=[1, 2]),
                        colormap=[[1, 124, 252, 0], [2, 0, 102, 0]], astype='u8')

Image(masked.export_image(bbox=area['extent'], size=[1200,450], f='image'))
<IPython.core.display.Image object>

Visualising results on map

m = gis.map('New Delhi, India', 10)
m.add_layer(masked)
m.legend = True
m

How-much-green-is-Delhi_output.PNG

Here, the pixels with light green color having a value of "1" represent agricultural land whereas the pixels with dark green color having a value of "2" represent forest area/ tree cover.

Area Derivation

Now we will calculate the total green cover of Delhi. It is important to note that pixels belonging to agricultural land as well as forest area/ tree cover are considered as green cover in the calculation.

pixx = (delhi_clip.extent['xmax'] - delhi_clip.extent['xmin']) / 1200
pixy = (delhi_clip.extent['ymax'] - delhi_clip.extent['ymin']) / 450

res = masked.compute_histograms(delhi_clip.extent,
                               pixel_size={'x':pixx, 'y':pixy})
numpix = 0
histogram = res['histograms'][0]['counts'][0:]
for i in histogram[1:]:
    numpix += i
sqmarea = numpix * pixx * pixy # in sq. m
acres = 0.00024711 * sqmarea   # in acres
HTML('<h3>Total Green cover in Delhi is ~ <i>{}%</i> of the total \
     geographical area of Delhi.</h3>'.format(int((acres/419004.17)*100)))

Total Green cover in Delhi is ~ 26% of the total geographical area of Delhi.

plt.title('Green Cover', y=-0.1)
plt.pie(histogram, labels=['Non green','Agricultural Land', 'Forest Cover']);
plt.axis('equal');
<Figure size 640x480 with 1 Axes>

The pie chart clearly shows the distribution of agricultural land, forest cover and non green areas like urban areas, water, etc. Here, the non green areas contribute around 74% of Delhi's land and remaining 26% goes with forest cover and agricultural land.

Conclusion

In this study, green cover of Delhi is calculated using Landsat 8 imagery for 21 October, 2022. Normalised Difference Vegetation Index (NDVI) is used for the calculation of green areas, which is a well known and widely accepted spectral index for vegetation studies.

NDVI is used here compared to the other spectral indices because it is easy to compute and requires only two bands. Pixels with NDVI values greater than 0.4 and less than 0.5 are considered as agricultural land, parks, etc and the pixels with more than or equal to 0.5 NDVI values are considered as tree cover/ forest areas.

Finally, the area of the pixels of agricultural land and forest class is calculated to estimate overall green cover of Delhi.

The study shows how the green land cover of an area could be easily computed in few lines of code using Esri's predefined NDVI layer and Landsat 8 imagery from ArcGIS server.

References

[1] https://fsi.nic.in/forest-report-2021-details

[2] https://economictimes.indiatimes.com/news/india/environmentalists-question-govt-report-on-delhis-forest-cover/articleshow/92164943.cms?from=mdr

[3] https://earthobservatory.nasa.gov/features/MeasuringVegetation

Your browser is no longer supported. Please upgrade your browser for the best experience. See our browser deprecation post for more details.