Last updated: 2022-10-24

Purpose

The purpose of this code is to adjust raster files for mapping - by reducing their resolution and clipping them to the state boundary. Note that this is simply for ease of creating figures. When using rasters in a downstream research application, it would be preferable to keep the resolution intact.

Load libraries

library(sf)
library(USAboundaries)
library(raster)

Load map files

# load one raster to get CRS info
pa_or <- raster("../output_big/rasters/test/CDL_2012_42.ld50_or_ha_bil.tif")

# PA border - download via USAboundaries
pa <- us_states(states = "Pennsylvania")
st_crs(pa) # check coordinate reference system
pa_GRS80 <- st_transform(pa, projection(pa_or)) # re-project to CRS of rasters
pa_poly <- as(st_geometry(pa_GRS80), "Spatial")

# download CDL raster
cdl <- raster('../data_big/cdl/CDL_2012_42/CDL_2012_42.tif') %>%
  mask(pa_GRS80)

# load pesticide rasters, reduce resolution, and clip to PA boundary
pa_or <- raster("../output_big/rasters/test/CDL_2012_42.ld50_or_ha_bil.tif") %>%
  raster::aggregate(fact = 2, fun = "mean") %>%
  mask(pa_GRS80)

pa_gly <- raster("../output_big/rasters/test/gly/CDL_2012_42.kg_ha.tif") %>%
  raster::aggregate(fact = 2, fun = "mean") %>%
  mask(pa_GRS80)

pa_mef <- raster("../output_big/rasters/test/mef/CDL_2012_42.kg_ha.tif") %>%
  raster::aggregate(fact = 2, fun = "mean") %>%
  mask(pa_GRS80)

pa_imi <- raster("../output_big/rasters/test/imi/CDL_2012_42.kg_ha.tif") %>%
  raster::aggregate(fact = 2, fun = "mean") %>%
  mask(pa_GRS80)

# survey coverage
pa_cov <- raster("../output_big/rasters/test/CDL_2012_42.coverage.tif") %>%
  raster::aggregate(fact = 2, fun = "mean") %>%
  mask(pa_GRS80)

# check that CRS for PA boundary and CDL are the same
pa_cdl <- cdl[[1]]
projection(pa_or) == projection(pa_cdl)
projection(pa_or) == projection(pa_GRS80)

Export map files

writeRaster(cdl, '../data_big/cdl/CDL_2012_42/CDL_2012_42_clip.tif')
writeRaster(pa_or, '../output_big/rasters/test/CDL_2012_42.ld50_or_ha_bil_clip.tif')
writeRaster(pa_gly, '../output_big/rasters/test/gly/CDL_2012_42.kg_ha_clip.tif')
writeRaster(pa_mef, '../output_big/rasters/test/mef/CDL_2012_42.kg_ha_clip.tif')
writeRaster(pa_imi, '../output_big/rasters/test/imi/CDL_2012_42.kg_ha_clip.tif')
writeRaster(pa_cov, '../output_big/rasters/test/CDL_2012_42.coverage_clip.tif')

This R Markdown site was created with workflowr