2

If we have a raster, say of integer elevation data for a country, and one polygon shapefile, say of 300 river basins in that country, with a unique name for each, how would we most easily get output like this for them all?

basinID, gridcellelev
a, 320
a, 321
a, 320
b, 17
b, 18
b, 19

The most burdensome way seems to be filtering/converting the single shapefile into 300 shapefiles, clipping the raster 300 times into 300 uniqueID rasters, reading them back in, generating individual tables for each basin, then combining them all together.

On the other hand the ideal way seems to be skipping the file generation, not saving the xy data, and creating the same table from just using one raster and one shapefile - by iteratively selecting the cells within a basin, stamping them with the basinID, creating a table, losing the coordinates, and keep repeating and appending that table until the 300th basin.

I'm not looking for any statistics, just the raw data listing of the grid cell elevations that would have been part of some standard clip method. I believe the attribute table that comes with the clipped raster from ArcMap is the counts/frequency of the cells. That output table format works for me too.

I don't know to minimally reproduce a raster and a polygon shapefile, so I'd just be grateful for any tips/libraries/functions/examples. If starting with R, here's a starting point:

library(tidyverse)
library(raster)
library(rgdal)
library(sf)

elev_raster <- raster("spain_elev_meters.tif") #integer raster
basins <- readOGR("spainbasins.shp", "spainbasins") %>% st_as_sf() #unique basin ID column: `basinID` 

I'd prefer to do everything in R, but happy to try any batch methods in ArcGIS as well (I have 10.6 with Spatial Analyst, but not Pro).

  • "I believe the standard raster clip output that comes out of ArcMap is not the raw data but the counts/frequency of the cells." No, you'll get a raster which is equal the original one but clipped on the area defined by one or more polygons. – umbe1987 Jul 19 at 14:29
  • Oh right, good point @umbe1987 - I'll revise - though the accompanying attribute table that I'm looking for seems to be the counts. – dbo Jul 19 at 14:31
  • The attribute table of a raster generally gives you only a list of unique pixel values witout repetition along with their respective count in that raster. – umbe1987 Jul 19 at 14:36
  • @umbe1987, either would be fine - the frequency as you mention or the raw list of each cell's elevation - as long as each row has the basinID. In R I'm familiar with how to quickly convert one to the other. – dbo Jul 19 at 14:40
4

While this is not a complete answer, as I do not have your exact files etc., it should serve as a close basis for what you are trying to do.

setwd("D:/blah/")
library(velox)
library(raster)
library(rgdal)
library(sf)
library(tidyverse)
elev_raster <- raster("ElevationRaster.tif")
basins <- readOGR("basin.shp", "basin") %>% st_as_sf() 
basins_ID <- data.frame(c(1:length(basins$geometry)),as.character(basins$name),stringsAsFactors = F) # create a frame to link geometry ID to the basin name
colnames(basins_ID) <- c("ID","Name") # a bit of house-keeping to make the frame easier to understand
vlx_data <- velox(elev_raster) # velox rasters are super fast for extracting. Normal raster-package ones are sloooow.
extracted_data <- vlx_data$extract(basins,df=T,small=T) # extracting with velox is easy
colnames(extracted_data) <- c("ID","Elev") # house keeping again
merged_data <- merge(extracted_data,basins_ID,by="ID") # and now we merge the data tables to make nice looking dataset
head(merged_data)

This results in this nice dataFrame:

  ID Elev      Name
1  1 31.050    a
2  1 14.550    a
3  1 19.840    a
4  1 63.226    a
5  1 86.181    a
6  1 77.979    a
  • awesome. thank you! very glad to know about velox now, too. – dbo Jul 19 at 15:52
2

If you're comfortable with arcpy, you can try this script.

Basically, it loop through each basin, clip the input raster and write in the output table all the values found in the clipped raster along with the name of the current basin.

import arcpy
from arcpy.sa import *

arcpy.env.overwriteOutput = True

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# CHANGE THESE TO FIT YOUR NAMES!!!
# COMPLETE PATH AND NAME OF YOUR INPUT RASTER
in_raster = r"C:\Umberto\corsi\umbe_arcgispro\GDEM30m_clip02.tif"
# COMPLETE PATH AND NAME OF YOUR BASIN SHAPE
basins = r"C:\Users\minorau\Documents\ArcGIS\Default.gdb\basins"
# NAME OF THE FIELD YOU WANT TO PRESERVE FROM THE BASIN ATTRIBUTE TABLE
basin_field = "name"
# OUTPUT PATH TO STORE YOU FINAL TABLE (I USED A GDB)
output_gdb = r"C:\Users\minorau\Documents\ArcGIS\Default.gdb"
out_table = arcpy.CreateTable_management(output_gdb, "output_table")

arcpy.AddField_management(out_table, basin_field, "TEXT")
arcpy.AddField_management(out_table, "Value", "LONG")

basins_lyr = arcpy.MakeFeatureLayer_management(basins, "basins_lyr")

# loop through each basin
with arcpy.da.SearchCursor(basins, [basin_field]) as b_cur:
    for basin in b_cur:
        arcpy.SelectLayerByAttribute_management(basins_lyr, "NEW_SELECTION", """{0} = '{1}'""".format(basin_field, basin[0]))
        raster_clip = ExtractByMask(in_raster, basins_lyr)
        ##Build attribute table for single band raster dataset
        ##Overwrite the existing attribute table file
        #arcpy.BuildRasterAttributeTable_management(raster_clip, "Overwrite")
        temp_result = arcpy.TableToTable_conversion(raster_clip, output_gdb, "output")
        with arcpy.da.SearchCursor(temp_result, ['Value']) as r_cur:
            for val in r_cur:
                with arcpy.da.InsertCursor(out_table, [basin_field, "Value"]) as i_cur:
                    i_cur.insertRow((basin[0], val[0]))

Here is a printscreen of the results and the inputs I've used to test it:

enter image description here

  • thank you! I'm most comfortable with R but I will learn from your example. – dbo Jul 19 at 15:51
  • You're welcome, happy you found a solution to your problem! – umbe1987 Jul 19 at 15:52
-2

How about "Convert Raster to Polygon", then do a "Spatial Join" of those new elevation polygons with your basin polygons. Do it so the centroids of the elevation polygons being inside of a basin polygon as the match method.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.