Workflow for updating raster using polygons of different attributes

Question

I require a specific output but it is proving difficult for me to work out a workflow.

I have:

1) polygons of new housing developments (356, non overlapping)

2) a 25m raster of land cover

I want to update my land cover raster (in areas where polygons overlay) to a new mix of land cover classes (predominantly built-up but with some improved grassland etc).

This mix is different for each polygon.

I have added to the polygon attribute table the number of cells of each land cover class that I want the polygon to represent when converted to raster (e.g. 18 built-up, 1 arable, 1 improved grassland, as the the polygon will be represented by 20 cells).

Furthermore, the cell values (land cover classes) need to be randomly allocated within the area covered by each polygon!

Any suggestions! Python or R code solutions welcome. I am more familiar with R but can run Python code if it doesn't need much adapting by me!

I am using ArcGIS 10.1 license type advanced. Windows 7 64-bit i7 processor

Answer 1

My solution was to convert the polygons to raster (snapped to the land cover map and with the same resolution) using a unique polygon ID field for raster value.

Then I converted this raster to points, and ended up with nearly 35,000 points, with attributes POINTID and GRID_CODE. POINTID is a unique cell identifier and GRID_CODE matches the original polygon ID. I joined the original polygon attribute table to the new points table using the GRID_CODE and polygon ID. This then gave me an attribute table like this for the points layer:

POINTID GRID_CODE cell_count cells_buil cells_arab cells_impr cells_brdl cells_semi
1       1        45          9          8          1          1          0          0
2       2        45          9          8          1          1          0          0
3       3        45          9          8          1          1          0          0
4       4        45          9          8          1          1          0          0
5       5        45          9          8          1          1          0          0
6       6        45          9          8          1          1          0          0

So I have the number of cells of each land cover class (fields 4-8) that I want to distribute among the points (cells) with the same GRID_CODE (e.g. in this case 45).

The rest I did in R.

I read this table into R and used the following code to create a land cover class field:

library(foreign)

points2<-read.dbf("C:\\temp\\housing_points.dbf")

# add landcover attribute field filled with value for builtup (16)

points2$lcover<-rep(16, nrow(points2))

# get list of polygon IDs without duplicate rows for looping through
id.list<-unique(points2$GRID_CODE)

#create empty dataframe
points.new<-points2[0,]

for (i in id.list){
  # select all points belonging to each polygon, one polygon at a time
  polygon.pnts<-points2[points2$GRID_CODE==i,]

  # extract the number of cells of each land cover class within that polygon
  cells.arab<-polygon.pnts[1,5]
  cells.impr<-polygon.pnts[1,6]
  cells.brdl<-polygon.pnts[1,7]
  cells.semi<-polygon.pnts[1,8]

  # create vector of new land cover classes (4 = arable, 5 = improved grassland, 2 = bleaf, 6 = semi-natural grassland)
  # will use to replace land cover later on
  new.lcover<-rep(4,cells.arab)
  new.lcover<-append(new.lcover, rep(5,cells.impr))
  new.lcover<-append(new.lcover, rep(2,cells.brdl))
  new.lcover<-append(new.lcover, rep(6,cells.semi))
    if (length(new.lcover) > 1) {
      new.lcover<-sample(new.lcover)
    } # shuffle the land cover values if more than one cell (point) being replaced

  # get total number of cells (points) of new land cover classes to know how many random rows are needed
  total.other<- sum(cells.arab, cells.impr, cells.brdl, cells.semi)

  # pick out relevant number of random rows, saving just the POINT_ID
  random.rows<-sample(polygon.pnts$POINTID, size = total.other)

  # overwrite built up land cover class for selected rows with other land cover classes (column 9)
  polygon.pnts[polygon.pnts$POINTID %in% random.rows,9]<- new.lcover
# join into one dataframe for exporting
  points.new<-rbind(points.new, polygon.pnts)
}

write.csv(points.new, "c:\\temp\\points_new.csv")

This took seconds to run in R

I then joined this csv file back to my points layer based on POINTID and converted to raster using the lcover field as the value.

Hope this helps someone in the future!