Raster analysis

Contents

• Part 1: Terrain Analysis
  • Objective:  Perform interpolation to create a digital elevation map (DEM), calculation of slope/aspect, and reclassification of rasters to locate the best place for a solar farm.
• Part 2: Zonal Operations
  • Objective: Perform zonal statistics.
• Part 3: Crime Density (Point Density)
  • Objective: :  Calculate the density of crimes in the city of Worcester.

---

Return to Geog 210 page

---

Part 1: Terrain analysis

In this part we want to find the best place for a solar farm.  The selected site must be on a relatively flat area (slope 0-5º of inclination) and where it gets the most solar radiation during the day (a south facing area).  In addition, it must be a single parcel larger than 60 acres. 

Note: Keep track of the procedures, you’ll have to make a cartographic model at the end.

Data: You don’t have direct access to a Digital Elevation Model (DEM) of your area of interest (AOI) but a good friend has provided you with an outline of the area, AOI (North of South Hadley) and a regular sampling of points with their elevation values (Elevation_Points).

First, we need to create a DEM from the point data through a process called Interpolation. 

Interpolation predicts values for cells in a raster from a limited number of sample data points. It can be used to predict unknown values for any geographic point data, such as elevation, rainfall, chemical concentrations, and noise levels.

– ArcGIS Help

We’ll use a method called Inverse Distance Weighted or IDW:

This method assumes that the variable being mapped decreases in influence with distance from its sampled location. E.g., when interpolating a surface of consumer purchasing power for a retail site analysis, the purchasing power of a more distant location will have less influence because people are more likely to shop closer to home

– ArcGIS Help

Task 0 – Preliminary – Activating Spatial Analyst extension

ArcGIS needs an extension called Spatial Analyst to work with rasters.  We have the extension installed but before using it we need to activate it.  If it’s not activated, you’ll get an error message when using some raster tools.

After launching ArcGIS Pro, before creating a map or opening an existing map, click Settings, on the left, then Licensing …  OR if you already loaded a map go to the Project tab (first tab on the left) \ Licensing. In the box “ArcGIS PRO extensions” verify that Spatial Analyst is licensed, If YES, you’re OK, keep working. If NOT, click the button “Configure your licensing options” below and check the checkbox next to Spatial Analyst.  Click OK and continue working.

Note: Spatial Analyst is a very old software that ESRI hasn’t updated in a long time so it cannot work with blank spaces in the paths or file names. Use “_” underscores if necessary.

Task 1 – Creating DEM and calculating Slope and Aspect

1. Create a new map project in a new working folder containing the required dataset.
2. Add the vector file Elevation_Points located in Terrain.gdb to your map and open the attribute table; notice the column containing elevation data in meters.
3. In the Analysis tab click the red Tools icon and in the Geoprocessing pane on the right, click on Toolboxes on top, then go to Spatial Analyst Tools\Interpolation\IDW

1. Enter Elevation_Points as input
2. Z value (this is the elevation value to be used in the DEM) should be set to the elevation field
3. Output raster, save it as DEM_AOI (no empty space in name)
4. Output cell size enter 30 (meters, don’t enter the units; units depend on the map units, meters)
5. Leave all the other fields as they are and click Run.

Study the resulting raster.  Verify its properties and make sure the cell size is 30 (Right click\Properties\ Source\Raster Information)

Question 1. How many columns does this DEM have?

The raster will show in the map window with about nine discrete classes.  Please note the color changes and the values in the legend carefully and study the terrain.  Can you see which areas are lower in elevation? Where are the higher elevations?

Question 2.  In general, in what direction is the elevation going up?
From North to South — South to North — East to West — West to East

Now we need to calculate the Slope.  This is the inclination of the terrain, or rate of change in elevation by change in distance.  You can express it as a percent change or in degrees of the angle θ (see figure).

• % Slope = (a/b) x 100   where “a” is the change in elevation (Rise) and “b” is the change in distance (Run)
• Slope in degrees = Arctan (a/b) (Arctan is the inverse tangent function of the ratio a/b; remember that a/b is Tan θ)
• Notice that a 100% slope is a 45-degree angle, larger angle will give you a slope > 100%

Knowing the elevation of each cell and their size (in our case 30m) the software calculates all the slopes (in either % or °) around each cell and outputs the steepest slope for the center cell. We’ll use slope in degrees.

4. Go to the Analysis tab click on the red Tools icon.

1. In the Geoprocessing pane click on Toolboxes on top then go to Spatial Analyst Tools\Surface\Slope
2. Enter DEM_AOI as input
3. Enter Slope_degrees as output
4. Output measurements as DEGREE
5. Leave Z factor at 1 (the vertical, height, values are in meter, same as the XY coordinates) and Run.

Study the result. You should get values between almost flat (0 degrees) to around 17 degrees.  If the software is showing you values up to ~90 degrees, go to the symbology pane of the layer and choose a Classify symbology with Natural Break 5 classes.  The legend should now show you the minimum and maximum values.

Now we need to find the areas that are facing the sun.  In the northern hemisphere, the sun is always on the southern side of the sky above us (remember that the sun shines perpendicular in tropical regions).  We need to find areas where the slopes are facing south.  This is the Aspect of the slope.

Using a compass as a guide, and knowing that a circle has 360 degrees around, we assume that north is 0º, east is 90º, south is 180º, and west is 270º (see figure).

5. Go to the Analysis tab click on the red Tools icon

1. In the Geoprocessing pane click on Toolboxes on top then go to Spatial Analyst Tools\Surface\Aspect
2. For input use the DEM_AOI (NOT the slope raster)
3. Save the output as Aspect_AOI.

Again, study the results: it should be very colorful. Notice the legend and copy the range of angles shown for South ________________ (copy these values). All aspects between those values are facing south and good for our solar panels.   These directions should be represented light blue in the CP, similar to the graph above.  Look at the map and try to see which areas have these colors. Compare it with the original DEM. Can you see the south facing areas in the DEM?

Notice also that there is a class for “flat” areas with a value of -1. Flat areas have no slope orientation.

Task 2 –Reclassification and overlay

Now we need to put those two themes (slope & aspect) together.  However, even though they both have units of degrees, they represent different things.  One represents degrees of inclination in the vertical direction; the other represents degrees of orientation from the north (in the horizontal direction).  We are going to separate the “Good” areas, just as we did with vectors in Exercise 6 land suitability.  This is called Standardization, i.e. taking layers that represent different things and making them compatible. We do this,  in raster analysis, with a function called Reclassification.

Reclassification assigns new numbers to our rasters according to our criteria.  In the case of the slope we are going to reclassify all the values with a slope of 0-5º  (i.e. flat) to have a value of “1”, all the other values will have a value of “0”.  For the aspect, we’re going to assign a value of “1” to our south facing slopes, and “0” to the rest.

6. Go to the Geoprocessing pane (Analysis tab\Tools), click on Toolboxes on top, then go to Spatial Analyst Tools\Reclass\Reclassify

1. Input: Slope_degrees (NOTE: use the SLOPE raster!)
2. Reclass field: Value
3. In the Reclassification table, make sure the first Start value shows the first class from the lowest value (0.008145, or just “0”) in the Start column
4. Enter a value of 5 in the End column of the first line, leave the New Value for this first class as 1
5. For the second class, enter 5 for start and 99 for End, change the New value to “0” (Any number higher than 17.32 will be good for the End field, 99 or 9999 is a “save” easy number)
6. Leave the NoData as it is
7. Left-click on the “extra” lines and press the delete key in the keyboard to remove them. Make sure there are only 3 rows: 0-5, 5-99, and NoData
8. Output raster Good_slopes. Click Run.

You should get a binary (0/1) map.  It’s important that you get a raster with only 0 and 1 value (this is called a Binary or Boolean raster).  If you got values 1 and 2, you didn’t assign the classes correctly, do it again.

Study the output; everything with a value of 1 is Good (relatively flat).

Since this is a raster with integer values, it does contain an attribute table, open it and check its content.

Hint: Good_slopes should have 29971 pixels with value of 1.

Now try to open the attribute table of the DEM, or the slope, or the aspect raster. You can’t, do you know why?

7. Now reclass the Aspect raster.

1. Input: Aspect_AOI (NOT DEM)
2. Reclass field: Value
3. In the Reclassification table enter, in the first row, a value -1 to 157.5, change the New Value for this first class to 0
4. For the second class, 157.5 to 202.5, change the New value to “1”
5. For the third class, 202.5-360, change the New value to 0. 
6. Leave the NoData as it is and delete the extra rows. You should end up with 4 rows
7. Output raster Good_aspects. Click Run.

Again, you should get a binary (0/1) map, if not, do it again. 

Hint: Good_aspects should have 5130 pixels with value of 1.

With vector layers, we would now do an intersection to find the sites where the two good conditions coincide.  In raster we use a trick from basic math:  all numbers multiplied by 0 become 0.  So we’ll just multiply our binary (0/1) maps and only in the places where cells with value 1 overlay another cell with value 1, there will be a 1 as a result.  See the little table below.  That is why it was important to have only 0’s and 1’s in our “Good” maps.

Good_Slopes		Good_aspects		Good_areas
0	x	0	=	0
0	x	1	=	0
1	x	0	=	0
1	x	1	=	1

8. Use the Raster Calculator (Geoprocessing pane \Spatial Analyst Tools\Map Algebra\Raster Calculator) to multiply Good_slopes and Good_aspects and name the output Good_Areas.

1. In the Raster Calculator pane double click on the 1^st layer (e.g. Good_aspects) then double click on the * symbol then double click on the 2^nd layer.  You “build” the following equation “good_aspects” * “good_slope“
2. Name output: Good_Areas

9. Examine Good_Areas. Only the few pixels with value of 1 are suitable for our purpose.

Question 3. What’s the total area of “Good_areas” in acres? 
Hint 1: Only pixels with value of 1 (suitable). 
Hint 2: Check the properties of the final raster to find out the area of one pixel, and the units of the map.
Hint 3: Knowing the area of 1 pixel, the attribute table tells you the number of pixels, so calculate the total area.  This area will be in map units square.   Convert those units to acres.

Very small patches of land are not good enough for a solar farm, so the total area calculated is not very realistic.  One final thing we need to do is to dis-aggregate all the separate patches (groups of touching pixels) with value 1 (equivalent, in the vector world, to converting a multipart polygon to single part).  Then we can extract the ones larger than certain critical size (in our case, how many acres? ______)

10. In the Geoprocessing pane launch \Spatial Analyst Tools\ Generalization\ Region Group

1. Input Good_Areas
2. Output Regions_good_areas
3. Number of neighbors FOUR (this decides which “touching” cells are counted as a group: FOUR – only cells touching up/down and left/right; EIGHT – same as before plus the diagonal cells)
4. Leave the other options with the default values and click Run.

Display Regions_good_area using a symbology of Unique Values and study its attribute table.  Notice the large number of new groups (regions or zones) created.  They have all been assigned a new Value, but they still have the information of suitable/unsuitable (1/0) in the column LINK, and, of course, the Count column tells you how many pixels there are in each region, therefore you can calculate the area of each patch of land.

HINT: there should be 328 regions

Question 4. How many groups of cells (regions) are SUITABLE? (Hint: 1’s)
Question 5. How many suitable parcels (regions) are larger than 60 acres?
(Hint: you need to calculate how many cells are in 60 acres)

In theory, we would also need to consider other factors like accessibility (distance to roads), environmental protection (buffers around water bodies), etc. but that’s a problem for another day.

Question 6. Make a cartographic model and upload to Moodle.
Hint: it’s very simple.
Hint 2: End it with the creation of Regions_good_areas

11. Save your map. It’ll be used in the next part.

---

Return to Top of page

Return to Geog 210 page

---

Part 2: Zonal operations

Task 1 – Zonal statistics.

After you explored the suitability of the area to have solar panels, a dispute among the neighbors has erupted.  They want to know how much of their land is suitable for solar panels and whether their homes are located in suitable areas.

To solve these problems, we’ll use a Zonal operator to evaluate Good_Areas using the neighbors’ parcels (polygons) and houses (points) as “zones”

1. Open the map from the previous part and add the layers AOI and Houses to your map (If you just have a new blank map, bring also Good_Areas and DEM_AOI)

The houses theme contains the information of the owner, but the parcels shown in AOI don’t have owner’s name (check both layers’ tables).  Therefore, we have to join these two.

2. Do a spatial join (Intersection) on AOI (target feature) with the houses theme (join feature).
3. Study the attribute table of AOI and make sure the names of the farmers are there.
4. Go to Geoprocessing pane \ Spatial Analyst Tools \Zonal \Zonal Statistics as Table (read the help of the tool to learn what it is that you’re doing)

1. Input raster or feature zone data (i.e. the file that defines the “zones”) – AOI
2. Zone field – Name (aoi_AddSpatialJoin.Name)
3. Input value raster (i.e. the raster to be evaluated) – Good_Areas
4. Output table – ZonalStats_Parcels
5. Statistics type – ALL

Study the output table.  It shows a series of statistics for all the pixels from Good_Areas that fell within the limits of each parcel.  Since the analyzed raster was a binary raster most of those columns don’t make sense, except SUM that shows you the sum of the values of all the pixels in the parcel. Since all the pixels only have values 0 and 1 their sum is the same as the total number of pixels with value 1. Notice that the column COUNT now includes both, pixels 0 & 1, it’s just a total pixel count.

Question 7. Which is the parcel (name) with the largest number of good pixels (with value 1)? 

Question 8. What is the total number of acres in the parcel above?
Hint 1: NOT just good pixels. 
Hint 2: don’t use the “AREA” column; the COUNT column tells you how many pixels there are, multiply these by the area of one pixel

Question 9. How many acres of good area are in this parcel? (Hint: <300) 

Finally, let’s redo this zonal statistic using DEM_AOI as our Input value raster.

5. Go to Geoprocessing pane \Spatial Analyst Tools\Zonal\Zonal Statistics as Table

1. Input raster or feature zone data (i.e. the file that defines the “zones” – AOI
2. Zone field – Name (aoi_AddSpatialJoin.Name)
3. Input value raster (i.e. the raster to be evaluated) – DEM_AOI
4. Output table – ZonalStats_DEM
5. Statistics type – ALL

Study the output table. Now, some of the calculated statistics do make sense.  You can now see for each parcel the minimum and maximum elevations (range is max minus min), the average elevation, and the standard elevation.  In this case the column SUM is the one that doesn’t make sense (can you tell why?).

Question 10. What is the average elevation for White parcel?

Question 11. What is the highest elevation found at Black parcel?

6. Save and close the map

---

Return to Top of page

Return to Geog 210 page

---

Part 3: Crime Density in Worcester (Point density)

Sometimes we want to know how many “events” (represented as points) happened within a particular area (like # houses/sq. miles). Rasters allow this calculation by counting the number of points that fall within a particular raster cell.

This part uses the same data in the GeoDB Worcester.gdb in exercise 6.  You can make a new map in \ex6 folder for simplicity.

1. Load Geocoded_Crimes (NOT the clipped version) and Worcester_City from the Worcester GeoDB into a blank map.
2. Go to Analysis\Tools\Geoprocessing pane \Spatial Analyst Tools\Density\Point Density

1. Input point features:  Geocoded_Crimes
2. Population field: leave as NONE
3. Output raster:  save as Crime_Density
4. Output cell size: 40  (the units of this distance is defined by the map units, meters, in this case)
5. Neighborhood: Circle.  Leave the Neighborhood settings as default
6. Area Units: Acres
7. Click Environments in the top of the dialog
   • Processing Extent: Same as Layer Worcester_City
   • Raster Analysis:  Mask Worcester_City  (this is necessary because without it, the software will compute density beyond the outline of the city)
8. Click run and wait for the density raster to display.

3. Represent the resulting layer using 10 categories and Geometrical Interval (this is a good example where an equal interval is a very poor representation). Use a dual tone color scheme.

Question 12. What is the highest crime density found in Worcester for this time period in crimes/acre?

---

Return to Top of page

Return to Geog 210 page