Remote Sensing using Python

I have been interacting with a potential oil-gas-energy client for whom I developed a proof of concept on using sound data at the oil rig sites. I have another meeting next week in which use cases around satellite imagery data will be discussed with the urban planning department of the firm. I want to learn the basics of GIS data analysis programmatically so that I am not exactly a neophyte by next week and come up with a business use case that’s relevant to the urban planning department. Thus, I have been in pursuit of a python library that’s powerful and strong enough to accomplish all the complicated tasks of remote sensing. That’s when I came across ArcGIS.

Oil and energy firms spend millions on the state of the art remote sensing software but most of it is used in oil exploration, measuring heat and moisture index, or hydrology studies.

Even though it is beyond my reach, I want to present a convincing factor that can govern a sustainable solution for the growth of urban areas in already choking cities. If the green cover is reducing at the expanse of human housing and other so-called development activities, then we need to look for alternative places and cities.

In this work, we will use the Python API of ArcGIS and explore the capabilities in finding the green cover of a city.

Let’s go!!

Problem Statement

Find out what percentage of London is green. Green areas include forests, parks, shrubs, and patches of green in the form of gardens that are visible from the satellites.

First things first, install the ArcGIS API. I would recommend that you do so in a virtual environment.

                     pip install arcgis

Check whether the install has finished successfully. The command below should execute in the Jupyter notebook without any error.

                        import arcgis

Create a free account at https://developers.arcgis.com/sign-up/ so you can use the userID and password for hitting the GIS API.

You can access the GitHub code at my repo here.

1. Imports

2. GIS connection

Connect to GIS and search for the multispectral Landsat layer.

3. Search for the UK

We will search for UK’s boundaries, it will have the boundaries at District, County, Census areas. My time in UK has taught me a thing or two, so I will go ahead with County level which is layer 1.

UK_items = gis.content.search('United Kingdom', item_type="Feature Layer", outside_org=True)[0]

UK_boundaries = UK_items.layers[1]

4. Find London

Find London’s using geocode and add it to Landsat’s extent object. London’s extent would be a key-value pair telling the width and height of the city of London on a map.

london_area = geocode("London, GB", out_sr=landsat.properties.spatialReference)[0]

landsat.extent = london_area['extent']

5. Find the longitude and latitude pair for London

In the UK_boundaries object, #7 belongs to London and we can extract all the pairs for latitude-longitude.

london = UK_boundaries.query(where='OBJECTID=7')
london_geom = london.features[0].geometry

london_geom['spatialReference'] = {'wkid':4326}

london.features[0].extent = london_area['extent']

6. Intersection of London in the Landsat data

selected = landsat.filter_by(where="(Category = 1) AND (cloudcover             <=0.2)", time=[datetime(2015, 1, 1), datetime(2019, 12, 31)],                          geometry=arcgis.geometry.filters.intersects(london_area['extent']))

Many datasets with such images would be present in Landsat, we need to filter the right one.

df_selected = selected.query(out_fields="AcquisitionDate, GroupName, CloudCover, DayOfYear",order_by_fields="AcquisitionDate").sdf

A dataframe will be released:

Let’s take objected = 2918532 from the above frame and see what does it do for us; feel free to grab any other id.

london_image = landsat.filter_by('OBJECTID=2918532')

apply(london_image, 'Natural Color with DRA')

7. NDVI(Normalised Difference Vegetation Index)

NDVI is an indicator of a plant’s health based on how a plant reflects different light waves. It does so by measuring the difference between near-infrared which plant strongly reflects and red light (which it absorbs). NDVI always ranges from -1 to +1.

Due to the surface reflectance property of the plants, NDVI metric can be used to identify the patches of green in a terrain.

ndvi_colorized = apply(london_image, 'NDVI Raw')

ndvi_colorized

8. Bringing everything together

Clip London’s latitude and longitude from the NDVI image above.

london_clip = clip(ndvi_colorized, london_geom)
london_clip.extent = london_area['extent']
london_clip

This is very neat! We have extracted London’s map from the Landsat image. Now we can find the green parts of the city.

9. The Green London

The clip of London’s map is the NDVI image, we can place a mask on it and find which parts are green.

10. Quantify the findings

Count the pixels in the image where the colour is green, one can also accomplish the same by drawing a histogram.

Total Green cover in London is ~ 39% of the total geographical area of London

If you have ever lived in London then you will know that it is very green and boasts of massive parks such as Regents, Hyde, Richmond, James, Victoria, Hampstead Heath, Battersea to name a few. The city is in the sweet spot where forest cover is just right(can be improved) for the number of people it caters(~10 million).

Conclusion

1. With python APIs, it is possible to do all that can be done with the software. Gone are the days when sophisticated tasks could only be performed using expensive software.
2. There is a learning curve associated with ArcGIS. The documentation will help you but you will have to spend a few hours to get around the concept.
3. London is a green city; the green cover should be closely monitored at repeated intervals to understand the dynamics of the city and to take corrective actions.

---

I like the parks of London and the free-roaming foxes the city has. The status quo should be maintained.

References:

API Reference for the ArcGIS API for Python – arcgis 1.8.3 documentation