This is the third entry in the satellite imagery analysis in Python series. For part I, click here, for part II, click here.

In this tutorial we will use land surface temperature (LST) as the data variable along with land cover information from the national (U.S.) database. The land cover information will allow us to create a relationship between land cover type and its respective heating (or cooling) contribution to the earth’s surface. Land cover is used in many applications ranging from algorithm development to military applications and crop surveying, not to mention applications in water management and drought awareness.

The Multi-Resolution Land Characteristics (MRLC) consortium produces the National Land Cover Database (NLCD) for the U.S. at 30-m resolution. The database defines 16 land cover classes for the lower 48 states (there are 4 more specifically for Alaska that we won’t use) that range from water, urbanization, forestation, and more. The NLCD from 2016 (the most recent version) can be downloaded from the following link:

https://s3-us-west-2.amazonaws.com/mrlc/NLCD_2016_Land_Cover_L48_20190424.zip

The file is quite large (1.4GB compressed and zipped). I recommend pulling the NLCD data into a GIS software and clipping the bounds to the particular area of interested (i.e. New York City in our case) - this will allow us to minimize the amount of processing needed between Python and the NLCD data.

In Python, the NLCD raster file needs to be read in using a geo toolbox. I chose to use rasterio because it is more familiar to me, however, other libraries are capable of handling the reading of the raster. The full methodology behind handling and plotting the NLCD raster file is quite involved and could be an entire tutorial, however, I will leave much of the explanation in the comments of the Python code. The breakdown for reading the NLCD is briefly explained below:

1. First, we read in the city boundary (NYC in this case) and set the bounding latitude and longitude points

2. Read the NLCD raster file to get its affine transform information - we will use this to inverse the NYC lat/lon boundaries to find the endpoints in raster file coordinates

3. Create a meshgrid for latitude and longitude in raster coordinates

4. Transform meshgrid to WGS84 coordinates using affine transform

5. Read NLCD data using lat/lon boundary indices

6. Read and delineate colors for NLCD bands

7. Map NLCD data

The full code for plotting NYC’s NLCD map is shown below, along with the corresponding figure output:

Above, we see that New York City is highly urban!

Below is a similar map of El Paso, TX. We can see the vast difference in land cover across the two cities. New York is highly urban, and while El Paso is urban as well, much of the outlying area is shrubbery and herbaceous desert-like land covers.

NLCD for El Paso, TX

Beneath El Paso, there is a region of darkness, which simply represents the change from the U.S. to Mexico. And since this is a U.S. land cover product, we do not have any data in that region. I’ve included a third NLCD map for the city of Phoenix, AZ below, which is fairly urban as well, with a lower population than NYC but great population than El Paso.

In the next section, we will dive into quantifying urbanization as well as other land covers to ultimately characterize temperature profiles that are inherent to particular cities and urban areas.

If we investigate the proportion of each land cover class for New York City (excluding 0 values and water), we can get an idea of the urbanization of land for the region. The following code snippet will print out the proportion of each land cover class:

The subsequent printout should be similar to the following:

-Developed, Open Space: 0.16
- Developed, Low Intensity: 0.19
- Developed, Medium Intensity: 0.29
- Developed High Intensity: 0.24
- Barren Land (Rock/Sand/Clay): 0.01
- Deciduous Forest: 0.04
- Evergreen Forest: 0.00
- Mixed Forest: 0.01
- Shrub/Scrub: 0.00
- Grassland/Herbaceous: 0.01
- Pasture/Hay: 0.00
- Cultivated Crops: 0.00
- Woody Wetlands: 0.02
- Emergent Herbaceous Wetlands: 0.02

These proportions add up to 1.0, which means that we have a proportion of land cover (minus water and NaN values) for the urban region. We can interpret this by using the Developed High Intensity land cover as an example. The high intensity proportion is 0.24, meaning that 24% of the land cover is highly urban for the region (NYC and some of New Jersey). Additionally, if we sum up all of the urban classes (Developed) we can state that the region is 88% urban - a striking statistic. The city of phoenix, by comparison, is only 45% urban. As we would expect, Phoenix is 47% shrub/scrub land - which is indicative of the desert-like land that surrounds it.

We can take this analysis even further by using the GOES-16 satellite pixel boundaries as clipping points to create a relationship between land cover proportion and land surface temperature (other variables as well, if desired). The method to follow is not necessarily a defined or exact process, but it can help us understand the influence that land cover has on heating of the earth’s surface.

First, we can start by bounding each GOES-16 pixel by its respective neighbors. This will allow us to clip the NLCD points by each satellite pixel’s rough bounding box. A snapshot of this happening in real-time is shown below, where a few clips have been computed and the bounding points are plotted within each GOES-16 pixel (where the pixel is defined at its center point):

Above, each satellite pixel center (black dot above) is bound by four red circles and will subsequently be filled-in by the colored rectangles. The colored rectangles represent the NLCD points contained within each GOES-16 bounded pixel (red circles). These boundings will allow us to calculate statistics relating GOES-16 data to ground parameters (NLCD). This is a powerful tool for analyzing the significance of land cover with, for our specific case, land surface temperature (LST). The code to replicate the plot above is given below, and the following plots and codes will demonstrate the implementation of the bounded colored rectangles.

If we wanted to map out how urban each satellite pixel was, we can use the bounding method shown above along with the standard satellite mapping technique, and end up with the following pixel mapping of urbanization:

In the image above, I’ve done a few things:

1. Technically, the figure is plotted with a shift to align the colors to the respective black dot satellite pixel centers. This is the true visualization of the satellite information. This is necessary because of the Python function ‘pcolormesh()’ which uses the top right corner as the input.

2. The urbanization percentage above uses only the Developed parameters of the NLCD

3. The urbanization is a linear combination of the developed parameters bound by each GOES-16 satellite pixel’s rough bounding

This method of combining the NLCD using the GOES-16 satellite pixel information will be the basis for the analysis that follow, which entails relating land cover proportion to LST data. This will give us an idea of how each land cover class affects surface temperature, indicating whether land cover is directly correlated to surface temperature or not.

The code for the urbanization of satellite pixels is given below:

Below is the historical weather from July 2019, which coincides with a heat wave in the north east. We will be using LST data from the heat wave period to correlate surface temperatures from different land cover classes that correspond to each NYC satellite pixel.

During the heat wave, we can investigate the heating of each land cover by comparing the proportion of land cover to the deviation from the city-averaged LST. This looks like the following for highly developed areas in New York City:

The figure above tells us that as we increase the high intensity developed areas of a satellite pixel, the surface temperature increases quite a bit from the rest of the land cover types. This has huge implications for urbanization and urban planning, indicating that if we continue to build high density cities, we could be producing warmer and warmer urban regions, which can have negative effects on the population of such cities.

Conversely, if we investigate less built-up land cover, we can see a slightly different result:

The low intensity urban land cover suggests that less built-up areas result in less net heating of the surface, something that could be useful when deciding how densely to build an urban area.

Furthermore, if we investigate the city of Phoenix, AZ and its highly-developed regions, the net heating of its urbanized surface can seen as an even more severe version of New York City:

The plot above suggests that Phoenix and its high intensity developed land cover results in a net 2K increase in LST.

If we look at the non-urban cropland area surrounding the dense regions of Phoenix, we can see the cooling effects of non-urban regions:

We can see that the effects of urbanization can result in an overall heating of the land surface, specifically when the area is considered high intensity. We can also the neutral effects of lower intensity urbanization and the net cooling effects of specific natural land cover (in the case of Phoenix: cultivated crops). These findings could be important for minimizing urbanization effects on global temperatures and beneficial for the general public, specifically around extreme heat events.

The code to reproduce the plots above is given below:

In this entry into the Satellite Imagery Analysis in Python series I explored land cover and its relationship to urban areas. Several cities and their land cover breakdown were introduced, along with a general method for clipping land cover into satellite pixel bins for comparison and analysis. Each city is different and unique because of its geography, natural landscape, and human influence; however, understanding the effects of urbanization can be beneficial for future urban planning and mitigating the negative effects of increased intensity of urban environments. In the last section of this tutorial, I introduced some results relating to the net heating of highly urban regions, something which could be further explored and used as a tool for analyzing the effects of cities on global temperature trends.

See More in Python and GIS: