ArcGIS Pro & Python — Analyzing Land Cover Change in Major Global Urban Centers
Populations around the world have increasingly shifted away from rural localities to urban centers. This trend has become so persistent that more than half of the global population now lives in urban places. This shift in geography demographics has the potential to cause unintended negative consequences, especially in regards to environmental issues.
Analyzing land cover change is a common way to observe the expansion of urban areas and the diminishment of vegetated areas to accompany population growth. Computing the change in land cover can provide regional planners, government authorities, and conservation organizations important data to develop cities and mitigate environmental hazards more effectively and sustainably.
Step One — Set the Working Environment
1–0. Geodatabases
Before downloading any data, I will create several geodatabases where I will save and store all data and subsequent layers that I will create. One geodatabase will be created for each of the five cities that I will analyze, and will be called:
“Phoenix.gdb”
“Singapore.gdb”
“Berlin.gdb”
“Buenos_Aires.gdb”
“Addis_Ababa.gdb”
To do this, I am going to open a new python notebook in ArcGIS Pro to run the python scripts and call it “Land_Cover_Notebook.ipynb.” The integrated python notebooks in ArcGIS Pro are based on Jupyter Notebooks, so they are intuitive and easy to use.
Here is my script to create a new geodatabase:
# Import system modules
import arcpy
# Create file geodatabase
arcpy.management.CreateFileGDB(r"C:\Users\cstile11\Desktop\Global Land Cover","Berlin.gdb")ArcGIS Pro also has the Python Window, which is a great option for running simple and quick scripts while seeing real-time results in the map above the window.
1–1. Coordinate Systems
After creating unique geodatabases, I have added a new map for each of the locations I am going to analyze. I now want to change the coordinate system of each map to match that of the location, and to make sure all of the coordinate systems are consistent with one another.
I am going to use the Universal Transverse Mercator (UTM) zones for my projected coordinate systems. UTM divides the world into 60 north-south zones, each 6 degrees of longitude wide. The zones are numberd consecutively, and coordinates are measured north and east in meters.
The ArcGIS Hub World UTM Grid is a useful interactive map that has all of the UTM zones of the world. I am going to use this web map to find all of my locations and which zone they lie within.
Here are the UTM Zones for each city:
Phoenix — UTM Zone 12
Singapore — UTM Zone 48
Berlin — UTM Zone 33
Buenos Aires — UTM Zone 21
Addis Ababa — UTM Zone 37
In the properties window of each map, I can find the UTM coordinate system that I want to use and change them. There are several UTM systems available in ArcGIS Pro; some are specific to continents and regions of the world. Since I am analyzing cities in various countries, I am going to use the “WGS 1984 UTM” system for each map. All locations, except for Buenos Aires, are in the northern hemisphere.
1–2. Land Cover Data
Land cover data is available for download from various websites — some global and some specific to countries. Instead of hunting for this data for each city online, I can navigate to the Living Atlas portal in ArcGIS Pro in the Catalog pane.
A quick search of “Land Cover” provides many results. Included in those results is “Global Land Cover 1992–2019” which is one raster layer that covers land cover for the entire globe. There is also a timeline to this dataset that allows me to select specific years and display those results, which is necessary to compute the changes between various years. I am going to add this data to each map.
This land cover raster layer is in 300 meter cell size resolution and 8 bit unsigned pixel type and depth.
Step Two — Edit Land Cover Data
2–0. Set the 1992 Timestamp
Since the land cover raster layer from the portal contains data from every year between 1992 and 2019, I need to timestamp each year to create a new layer so I can compare changes later.
To timestamp, I am going to create a new definition query that will display the data for the earliest time that is available, which is January 1, 1992.
Now the 1992 data is displayed in my map where I will export that data as its own layer.
2–1. Export the 1992 Raster
To export the 1992 data, I will use the “Export Raster” tool. Before running the tool, I will make sure that my map’s current display is at a 1:400,000 scale, which will be the same for every city’s map so I have consistency across all maps and layers.
This tool will also allow me to set the correct coordinate system as well as to clip down the layer to a more manageable size. An optional Python script, which I will show later, can accomplish the same task.
The new layer created has been called “Berlin_1992.gdb.”
2–2. Set the 2019 Timestamp
Setting the 2019 timestamp will be the same method as the 1992 data. I need to make sure this definition query is executed on the original raster layer that contains all years, not the new layer I just created for 1992.
2–3. Export the 2019 Raster
Just like with the 2019 data, I will use the “Export Raster” tool to make a new raster layer for the 2019 data. This layer will be called “Berlin_2019.gdb.”
I should now have two new rasters, one for 1992 and one for 2019. This process will be repeated for each city. These two rasters will be used for the analyses to detect changes throughout the years.
2–4. (ALTERNATE OPTION) Use Python to Clip & Change Coordinate System
Instead of setting the clipping extent and coordinate system while exporting the raster, here is an alternate option to accomplish those same tasks using a Python script.
The python script will use a “with” command and run the coordinate system first, then it will be able to execute the “Clip” tool. This time, I ran the script in the Python Window.
Here is my script to set the projected coordinate system and clip the raster for Berlin:
# Clip raster with correct coordinate system
with arcpy.EnvManager(outputCoordinateSystem='PROJCS["WGS_1984_UTM_Zone_33N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]', extent='320748.742880528 5775837.65853288 474065.083709664 5867894.41339754 PROJCS["WGS_1984_UTM_Zone_33N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",15.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'):
arcpy.management.Clip("Global Land Cover 1992-2019", "320748.742880528 5775837.65853288 474065.083709664 5867894.41339754", r"C:\Users\colin.stiles\Desktop\GIS Training\Tutorials\Global Land Cover\Berlin.gdb\Berlin_Land_Cover_1992_2019", None, "255", "NONE", "NO_MAINTAIN_EXTENT")Step Three — Computing Land Cover Change
With all land cover raster layers edited and ready to go for each city, I can now analyze the change in land cover between 1992 and 2019. There are several ways within ArcGIS Pro to analyze the change in rasters from two different years.
3–0. Alternate Tools
A common option that is used is the “Change Detection Wizard” which is a simple and intuitive tool that allows the user to preview results while running the wizard. The Wizard is a guided workflow that encompasses three steps (for categorical change): Configure, Class Configuration, and Output Generation.
There is also the “Compute Change Raster” geoprocessing tool, which accomplishes the exact same thing.
3–1. Executing the Python Script
Instead of using any of the GUI (graphics user interface) tools, I ran my analyses from the following Python script:
# Check Image Analyst extension license
arcpy.CheckOutExtension("ImageAnalyst")
# Compute land cover change from 1992 to 2019
out_raster_dataset = arcpy.ia.ComputeChangeRaster("Berlin_1992", "Berlin_2019", "CATEGORICAL_DIFFERENCE", "'Rainfed Cropland';'Herbaceous Cropland';'Tree or Shrub Cropland';'Mostly Cropland in a Mosaic with Natural Vegetation';'Mostly Natural Vegetation in a Mosaic with Cropland';'Closed to Open Canopy Broadleaved Deciduous Tree Cover';'Closed to Open Canopy Needleleaved Evergreen Tree Cover';'Closed to Open Canopy Needleleaved Deciduous Tree Cover';'Mixed Tree Cover';'Mostly Trees and Shrubs in a Mosaic with Herbaceous Cover';'Mostly Herbaceous Cover in a Mosaic with Trees and Shrubs';Grassland;'Sparse Vegetation';'Flooded Shrub or Herbaceous Cover';'Urban Areas';'Consolidated Bare Areas';'Bodies of Water'", "'Urban Areas'", "CHANGED_PIXELS_ONLY", "AVERAGE", "ClassName", "ClassName")
# Save the output
out_raster_dataset.save(r"C:\Users\cstile11\Desktop\Global Land Cover\Berlin.gdb\Berlin_Change")Checking the Image Analyst extension is not required, but having the extension license is required in order to run the ArcPy tool. For this Berlin example, the output raster was saved as “Berlin_Change.gdb” in the corresponding geodatabase.
Here is the output raster for Berlin:
By clicking on one of these cells, I can see the transformation that took place for the single cell, which is 300 square meters in size. For example, choosing a pink cell will tell me that the area went from Rainfed Cropland in 1992 to Urban Area in 2019. All of the land cover transformations are listed in the Contents pane next to the map.
Step Four — Analyze Land Cover Change Results with a Bar Graph
Creating a bar chart of the results is another useful way to better visualize the resulting changes in land cover. With the “Create Chart” tool in ArcGIS Pro, I can quickly create several types of charts visualizing my analyses.
4–0. Clean Attribute Data
Before creating the chart, I need to organize and clean the data. After opening up the attribute table on the raster and sorting by ‘Area’ in descending order, I can see that there is one row named ‘Other’ and one named ‘No Change’ containing all the pixels that did not experience any change between the years. I do not need these data for the chart, so I will delete them.
4–1. Create the Bar Chart
With the “Create Chart” tool, I am going to create a new bar chart. I want my 1992 land cover classes on the x-axis and the total area that was changed on the y-axis.
4–2. Edit the Chart Properties
After editing the chart series and formatting, here is the final bar chart that displays the type and amount of land cover from 1992 that changed to urban land cover in 2019:
Step Five — Compare All Cities
Each city ended up with different results, but the general theme is the same — in the years from 1992 to 2019, all five cities experienced significant urban growth. Much of this urban growth came at the cost of cropland, forested areas, and shrublands.
5–0. Phoenix Results
5–1. Addis Ababa Results
5–2. Buenos Aires Results
5–3. Singapore Results
5–4. Berlin Results
Conclusion
Using tools in ArcGIS Pro and Python, cities and regions have the ability to analyze changes in land cover throughout many years. As we have seen here, no two cities will have the same results. To take this analysis even further, we could compare year by year to identify which years experienced the most change, and why. These analyses can help cities plan for the future to grow sustainably and more efficiently, taking any guesswork out of development and prioritizing the preservation of natural landscapes.
Resources
ArcGIS Pro 3.0.2. Copyright 1995–2022. Esri, Inc.
