In this analysis, we use R to calculate zonal statistics of a raster. The data come from the 2006 National Land Cover Database.

First we relcassify the NLCD data into two classes, based on the value of the raster. In this case, I want to classify the pixels as to whether they are developed or undeveloped.

We then use a census tract polygon layer to calculate the proportion of each tract’s land area that is developed vs undeveloped.

Raster read in

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.3-2, (SVN revision 755)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/ozd504/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/ozd504/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
library(raster)

nlcd<-raster("~/Google Drive/classes/dem7093/GIS_class_2018/data/bexar_nlcddata.tif")
plot(nlcd, main="Original Raster image")

nlcd<-na.omit(nlcd)

Raster reclassification:

m<-c(1, 12, 0,13, 24, 1, 25, 95, 0)
rclmat<-matrix(m, ncol=3, byrow = T)
nlcdreclass<-reclassify(nlcd, rcl = rclmat)
plot(nlcdreclass$layer, main="Reclassified Raster image")

Load the census tract polygons

library(tigris)
## To enable 
## caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
## 
## Attaching package: 'tigris'
## The following object is masked from 'package:graphics':
## 
##     plot
satract<-tracts(state="48", county="029", year=2010)

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
satract<-st_as_sf(satract)
satract<-st_transform(satract, crs=  102009)

satract<-as(satract, "Spatial")
summary(satract)

Calculate zonal statistics

Now we find the mean of the binary reclassified data, to estimate the percent of each tract’s area that is developed.

First we reproject the raster to the same coordinate system as the census data. I get the projection string from ^ above.

nlcdreclass<-projectRaster(nlcdreclass, crs="+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
satract$pct_developed<-as.numeric(extract(nlcdreclass, satract, fun=mean)) #this takes a minute or two

satract<-st_as_sf(satract)

library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
satract<-satract%>%
  filter(complete.cases(pct_developed))

tractlines<-st_boundary(satract)

library(classInt)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`
m1<-satract%>%
  mutate(devcut=cut(satract$pct_developed,breaks=data.frame(classIntervals(var=satract$pct_developed, n=5, style="jenks")[2])[,1], include.lowest = T))%>%
  ggplot(aes( fill=devcut))+geom_sf()+
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
  geom_sf(data=tractlines,fill=NA, color="black")+
  guides(fill=guide_legend(title="Percent Developed"))+ggtitle(label="Percent of Land Area Developed, 2006", subtitle ="NLCD 2006")+theme(axis.text.x = element_blank(), axis.text.y = element_blank())+theme( legend.text = element_text(size = 14), plot.title = element_text(size = 18), strip.text.x = element_text(size = 16))

m1

I used this method in a paper I wrote in 2013, although in that paper I did the diversity of the land cover. That could be done like:

nlcdproj<-projectRaster(nlcd, crs="+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0")
satract$var_land<-as.numeric(extract(nlcdproj, satract, fun=sd)) #this takes a minute or two

m2<-satract%>%
  mutate(devcut=cut(satract$var_land,breaks=data.frame(classIntervals(var=satract$var_land, n=5, style="jenks")[2])[,1], include.lowest = T))%>%
  ggplot(aes( fill=devcut))+geom_sf()+
  scale_fill_brewer(palette = "Blues") + 
  scale_color_brewer(palette = "Blues")+
  geom_sf(data=tractlines,fill=NA, color="black")+
  guides(fill=guide_legend(title="Land Use Diversity"))+ggtitle(label="Standard Deviation - Land Use, 2006", subtitle ="NLCD 2006")+theme(axis.text.x = element_blank(), axis.text.y = element_blank())+theme( legend.text = element_text(size = 14), plot.title = element_text(size = 18), strip.text.x = element_text(size = 16))

m2

Exporting images for posters

I recommend using the ggsave() function to export images for use in other documents:

ggsave(m1, path = "~/Google Drive/classes/dem7093/GIS_class_2018/data/", filename = "sa_pct_devel.png", units = "in", width=8, height = 10)

And you can insert that image into a doc or presentation.