We will be quantifying land cover changes that occurred between 2004 and 2018 today, using vector cover type maps like the ones we worked with the first week of class, as well as their raster-ized version. We will be able to see how change detection is done using vector or raster data this way, with the same set of cover types.
These cover type maps are maintained by the San Diego Association of Governments (SANDAG), and they are primarily used for land use planning purposes. They are thus much more detailed in the developed land cover types than they are in their undeveloped cover types - most vegetated areas are just called "Landscape open space", and there are no categories for different habitat types. But, we will be able to tell when undeveloped land in 2004 has been developed, and what it has been converted to, which is sufficient for tracking conversion of undeveloped land to developed land. To keep the analysis relatively simple I condensed the 60+ different cover types in the original data into a set of 11, most of which are different types of human development.
We could use any two years for this exercise, but these two were chosen because we have NAIP images taken near the same time that the cover type maps were produced, and we will be able to make use of the NAIP images to assess the accuracy of the maps in our next exercise.
Let's start with the raster versions of the maps. Start up ArcMap and load lu_2004_rast and lu_2018_rast, which are the vector cover type maps converted to a raster format. I used 30 x 30 m pixels to be consistent with the size we have been using for our LandSat bands.
If you right-click on the raster maps and open their attribute tables, you'll see that the attribute tables have both a numeric code and a text label in them, which makes them easy to interpret. To help you keep track of which map is which, the label for land use is called "LU_2004" for the 2004 data, and "LU_2018" for the 2018 data. All we need to do is to compare pixel to pixel between the two maps, and count up how often the land uses are the same between the years, and how often they are different - this is called "cross tabulating".
Once we are done we will have a "transition matrix" that will look something like this (with only four of the cover types shown to keep the example simple - these numbers are not the ones you will get from the actual map, and are just for illustrative purposes):
| Cover in 2018 | ||||||
| Cover in 2004 | Residential | Landscape open space | Water body | Agriculture | Total for 2004 | |
| Residential | 900000 | 0 | 0 | 0 | 900000 | |
| Landscape open space | 1800 | 2700000 | 9000 | 18000 | 2728800 | |
| Water body | 1800 | 3600 | 450000 | 900 | 456300 | |
| Agriculture | 0 | 0 | 0 | 81000 | 81000 | |
| Total for 2018 | 903600 | 2703600 | 459000 | 99900 | 4166100 |
|
The cover type of a pixel in 2004 is indicated in the row labels, and the cover type of the same pixel in 2018 is indicated in the column labels. The data are expressed as square meters. To interpret the table, pick a cover type from 2004 in one of the rows, and read across the columns to see what the land cover has become in 2018. You'll see that all of the 900,000 square meters of residential area in 2004 was still there in 2018, and none of it became landscape open space, water bodies, or agriculture. There was a total of 2,728,800 square meters of landscape open space in 2004, of which 2,700,000 were still landscape open space in 2018, but 1800 square meters was converted to residential, and 18,000 square meters was converted to agriculture, and 9,000 square meters became water bodies (presumably due to habitat restoration, or changes in reservoir level). Similarly, water body area was still water body in 2018 (450,000 of 456,300 m2), but some was converted to residential (1,800 m2), some was converted to agriculture (900 m2), and some was converted to landscape open space (3,600 m2, probably due to ecological succession, or low water level in reservoirs). All of the agriculture from 2004 was still present in 2018.
We will calculate a table like this, using all 11 of the cover types in the cover type maps. ArcMap will give us the cross-tabulation, and we can add totals in Excel.
1. Cross tabulate the maps.
First we will need to generate the cross-tabulation of the raster maps in ArcMap.
Open the ArcToolbox, and navigate to "Spatial Analyst Tools" → "Zonal" → "Tabulate Area" - double-click to start the tool.
Select lu_2004_rast as the "Input raster or feature zone data", and select "LU_2004" as the "Zone field".
Specify lu_2018_rast as the "Input raster or feature class data", and "LU_2018" as the "Class field".
Put the output table in your monitoring.mdb, and call it "trans_rast".
The processing cell size will come from the lu_2004_rast layer (the first one you selected) - the tool will read the size of the pixel and use that to convert the counts of pixels into areas. The pixels are 30 x 30 m, but the units of measure used by the layer is feet, so the numbers we get in the table will be in square feet. How would you know this if I didn't tell you? See below...
Click "OK" to run.
When the analysis is complete you'll have trans_rast listed as a table in your Table of Contents. Note that it also sets the Table of Contents to "List by Source" mode, which shows the location of the files in your TOC, complete with drive letters and paths. You can't drag files around to re-order them in this mode - if you need to re-arrange files, you would need to switch back to "List by Drawing Order" mode, which you select by clicking on the left-most icon just beneath the "Table Of Contents" window title bar. Keep the TOC in "List by Source" mode for now, though, so we can export the table to an Excel file.
**How can you know the units used for trans_rast? The units are read from lu_2004_rast, and you can see what they are by right-clicking on that layer, opening its "Properties", and switching to "Source" - if you scroll down to the "Raster Information" section it gives the size of the cells, and the "XY Coordinate System" section lists the "Linear unit" as being "Foot_US" - so, each pixel is 98.4252 feet x 98.4252 feet, which is the same as 30 m x 30 m. Why, you ask, did I use a spatial reference that uses feet when I wanted to define pixels in nice, round metric units? The answer is, I wasn't thinking straight. But, having made the mistake it seemed a good opportunity to show you how to look up the sizes of pixels in a raster map, so teachable moment and all that.**
2. Export the file, open it in Excel. If you right-click and open the trans_rast table you just created (double-clicking will not work!), you will see it has a column for "LU_2004", which gives the cover type in 2004. The columns to the right are labeled "AGRICULTURE", "BEACH", etc., which are the category names for the 2018 file - since they were converted to column names, which cannot have spaces, the spaces are replaced with underscores, but otherwise they will match the categories from the 2004 file, and will be in the same alphabetical order.
We want to get these data into Excel so we can do some data summary, like so:
3. Calculate totals, convert to hectares, and calculate percents of rows. First we'll calculate totals for 2004.
Next we'll calculate the 2018 totals:
Now, we want to convert the numbers from square feet to hectares - the numbers are huge as square feet and difficult to interpret (they are also not metric units, which is not very scientific). There are 107639 square feet in a hectare, so if we divide the numbers in square feet by 107639 we will get area in hectares.
Now we will make one last version of the table that expresses the data as percentages of the row totals - since the row totals are the total area of each cover type in 1990, the percentage of row total will tell us what percentage of each land cover type in 2004 became each of the cover types in 2018.
If all went well, you'll see that 72.6% of agriculture in 2004 was still agriculture in 2018, and 45.7% of beach
was still in that cover type in 2018.
5. Check for tell-tale signs of error. There are two basic reasons for a transition from one cover type to another in this table: 1) a conversion from one type to another in the real world, or 2) an error in the 2004 and/or 2018 cover type maps. Some types of errors won't be apparent (e.g. a change from water body to landscape open space could be due to ecological succession, or it could be a mistake), but some types of transitions are so unlikely they are almost certainly an error. For example, residential or commercial areas in 2004 are very unlikely to become landscape open space in 2018.
6. Save the data. Save your Excel spreadsheet - call it "transition". This is the transition matrix for the land use changes that you'll use in your report.
Now we will see how to use vector polygons to do change detection work. Back to ArcMap...
1. Load landuse_2004_simple.shp and landuse_2018_simple.shp. These are land cover vector maps, like the ones we used in the first lab, except that they come from different years and I have simplified their cover type categories.
2. Do a spatial overlay (union) of the two shapefiles. The union operation overlays the features in both of the shapefiles, makes new polygons from the overlapping polygon outlines, and makes an attribute table for the resulting shapefile that has all of the columns from both of the attribute tables from the original files. Because of this you will end up with an attribute table that tells you what the land use was in 2004, and what it is in 2018. The area of each new polygon is calculated as well, so we can get a transition matrix by cross-tabulating the cover types, and adding the areas together.
To do the overlay:
From the ArcMap menu, select "Geoprocessing" → "Union".
Drop down the "Input features" box, and first select landuse_2004_simple, then do it again and select landuse_2018_simple.
Set the "Output Feature Class" to output to monitoring.mdb, and call the new layer "landuse_union".
Keep all other settings at default values, and click "OK" to run it.
This might take a few minutes - as long as the "Union..." progress indicator is running at the bottom of the map window be patient, it's still working.
Once it's done, you can zoom in to one of the areas with lots of lines, and turn the landuse_union layer on and off to see how it compares with the other two layers. You should see that all of the lines from both landuse_2004_simple and landuse_2018_simple are combined in landuse_union.
3. Export the attribute table. Right-click and open landuse_union's attribute table. You'll see that there is an FID for each of the layers you overlaid, and a column identifying the land use in 2004 and another for 2018 as well. If you scroll to the right you'll see there is also a Shape_Length and Shape_Area column, which give the perimeter and surface area for each polygon, respectively. The units are taken from the data set - since the units used in these layers is feet, the perimeter is in linear feet and area is in square feet.
The units are a bad choice, because the numbers are huge and hard to compare. We should also pick something metric, because this is Science! To get units in hectares, do the following:
Drop down the "Table Options" menu (upper left corner of the Table window), and select "Add field".
In the "Add Field" window, use "Area_ha" as the name, and select the "Type" as "Double". We will be calculating an area in hectares (ha), and we want to use double-precision floating point numbers so that small polygons don't get assigned an area of 0.
Click "OK" to create the new field (column in the table). You will now have a blank column, with <Null> appearing in all of the rows.
Right-click on the column label, Area_ha, and select "Calculate geometry". You will be warned that you are going to do something that can't be undone, which is fine - say "Yes". In the "Calculate Geometry" window, leave the property set to Area, but change the Units to hectares. click "OK". You will again be warned that you can't undo the calculation, which is still fine - click "Yes".
You should now have area in hectares.
Now that you have the area in hectares column we can export the table to an Excel file.
4. Open the exported table in Excel, and make the transition matrix. Switch back to Excel, and open the transition.xlsx file you made earlier (if it is not already open). We are now going to get a copy of the data from landuse_union.xls into this file:
You can close landuse_union.xls now, you have a copy of the data in transition.xlsx.
Now to make the transition matrix for the vector data we will use a PivotTable. If you haven't used PivotTables in Excel you're in for a treat - they make your life better.
You'll see that this creates a table in which the land use in 2004 is in the rows, the land use in 2018 is in the columns, and the total area of land in each combination of 2004 and 2018 land uses is in the body of the table. This is your transition matrix - the row and column labels sort automatically, and the units are already in hectares.
The tab is called "Sheet2", which isn't very descriptive - double-click on it and change it to "Vector", since this is the transition matrix for the vector maps.
The table you have is using Area_ha, so the units are already hectares and there's no need to do a conversion. Pivot Tables automatically give row and column totals, so there's no need to do those calculations either.
There are a small number of blanks for one year or the other - small differences in where the coast was drawn are responsible for this, because there is a little bit of area from 2004 that doesn't have any land use polygons to match it in 2018, and vice versa. We can hide these by dropping down the menus next to "Column Labels" and "Row Labels" and un-checking (blank)
To get the data expressed as a percent of row total we just need to change the display settings for the Area_ha field:
You can go back and forth between hectares and percent of row totals as needed. To view both at once you can copy the pivot table and paste it to cells below the current one:
You'll see that the values are similar to, but not identical to, the raster map version of these data. Converting between raster and vector maps is not perfect, and some amount of change in area can occur (how much depends on the size of pixels used - the distortion is greater the bigger the pixel size).
Save your transition matrix Excel spreadsheet, and a map file for today.