Raster Quiz
Josh Carr
2/20/2023
Set-up
library(raster)
library(tidyverse)
library(rgdal)
library(sf)
library(spdplyr)
library(maptools)
library(terra)Q1: Download the CropScape raster for Cache Valley in 2016. Replace
all pixels with a value of 0 (zero) with the value
NA, then compute the frequency table for this raster. What
was the most common land use category (after NA) in the
raster?
Code to load raster data.
cache <- raster("C:/Users/joshc/Documents/USU/2023/GEOG 4870/Module 6_Rasters/Quiz 6/CDL_2016_49005/CDL_2016_49005.tif")
spplot(cache)## Warning in lapply(xd, as.numeric): NAs introduced by coercion
The most common land use category was 141; Deciduous Forest.
Code to replace pixels with 0 value to NA.
cache[cache == 0] <- NADouble check 0 was replaced with NA by viewing minimum value.
minValue(cache)## [1] 1Hint: Look at the raster metadata here and review this week’s tutorial for reminders on how to download the data.
Code to load the frequency table.
freq(cache)## value count
## [1,] 1 41082
## [2,] 4 691
## [3,] 12 30
## [4,] 21 45817
## [5,] 23 3692
## [6,] 24 92755
## [7,] 28 1652
## [8,] 33 34813
## [9,] 36 355800
## [10,] 37 77071
## [11,] 43 1281
## [12,] 49 16
## [13,] 57 3
## [14,] 59 239
## [15,] 61 46904
## [16,] 66 56
## [17,] 67 63
## [18,] 68 10
## [19,] 70 1
## [20,] 111 22427
## [21,] 121 66079
## [22,] 122 52769
## [23,] 123 21041
## [24,] 124 7690
## [25,] 131 5098
## [26,] 141 897396
## [27,] 142 578479
## [28,] 143 7924
## [29,] 152 865163
## [30,] 176 90818
## [31,] 190 13812
## [32,] 195 43205
## [33,] 205 1850
## [34,] 229 15
## [35,] NA 1283028Q2: What percent of non-NA pixels have a value of 36
(alfalfa)? (Emphasis on the non-NA!)
percent_alfalfa <- length(cache[cache == 36]) / length(cache)
percent_alfalfa## [1] 0.07637209Alfalfa makes up 7.6% of non-NA pixels.
Q3: Crop the Cache Valley land use dataset to the
boundaries of the Logan municipality (not the extent
of the Logan municipality, but the actual shapefile
boundaries). You’ll need to extract the Logan shapefile from
the larger Municipalities shapefile found on the Utah GIS Repository
(HINT: spdplyr). Go to this
link and scroll down until you find the
Municipal Boudnaries:Shapefile download link. Upload your
plot on Canvas (2)
Upload Logan Munipality Shapefile.
lm <- readOGR(dsn = "C:/Users/joshc/Documents/USU/2023/GEOG 4870/Module 6_Rasters/Utah_Municipal_Boundaries/Municipalities.shp")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\joshc\Documents\USU\2023\GEOG 4870\Module 6_Rasters\Utah_Municipal_Boundaries\Municipalities.shp", layer: "Municipalities"
## with 254 features
## It has 18 fieldslm## class : SpatialPolygonsDataFrame
## features : 254
## extent : -114.044, -109.3116, 36.99999, 42.00117 (xmin, xmax, ymin, ymax)
## variables : 18
## # A tibble: 254 × 18
## OBJECTID COUNTYNBR NAME COUNTYSEAT SHORTDESC UPDATED FIPS ENTITYNBR
## <int> <chr> <chr> <int> <chr> <chr> <chr> <dbl>
## 1 1 12 Nephi 1 NEPHI CITY 2022/1… 54220 3040
## 2 2 27 Rockville 0 ROCKVILLE <NA> 64560 3075
## 3 3 17 Woodruff 0 WOODRUFF 2017/1… 85260 3040
## 4 4 21 Sigurd 0 SIGURD TOWN 2023/0… 68650 3110
## 5 5 08 Cleveland 0 CLEVELAND 2023/0… 14070 3030
## 6 6 29 South Ogden 0 SOUTH OGDEN 2023/0… 70960 3100
## 7 7 23 Wendover 0 WENDOVER 2016/0… 82730 3070
## 8 8 03 Richmond 0 RICHMOND 2022/1… 63680 3150
## 9 9 29 Riverdale 0 RIVERDALE 2022/0… 64010 3080
## 10 10 11 Brian Head 0 BRIAN HEAD… 2010/1… 08020 3010
## # ℹ 244 more rows
## # ℹ 10 more variables: SALESTAXID <chr>, IMSCOLOR <int>, MINNAME <chr>,
## # POPLASTCEN <int>, POPLASTEST <int>, UGRCODE <chr>, GNIS <chr>,
## # GlobalID <chr>, SHAPE_Leng <dbl>, SHAPE_Area <dbl>NL.sh <- filter(lm, lm@data$FIPS == "54990" )
L.sh <- filter(lm, lm@data$FIPS == "45860")
# North Logan:
extent(NL.sh)## class : Extent
## xmin : -111.8477
## xmax : -111.7708
## ymin : 41.75839
## ymax : 41.79492NL.extent <- extent(c(-111.8477, -111.7708, 41.75839, 41.79492))
NL.matrix <- matrix(st_bbox(NL.extent), nrow = 2)
# Logan:
extent(L.sh)## class : Extent
## xmin : -111.902
## xmax : -111.7802
## ymin : 41.67539
## ymax : 41.79736L.extent <- extent(c(-111.902, -111.7802, 41.67539, 41.79736))
L.matrix <- matrix(st_bbox(L.extent), nrow = 2)Reprojecion:
# Cache reprojection
crs <- '+proj=utm +zone=11 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'
cache <- projectRaster(cache, res = res(cache), crs = crs, method = "ngb") # takes a sec
# NL Transform
NL.sh <- spTransform(NL.sh, crs)
# L Transform:
L.sh <- spTransform(L.sh, crs)
cache@crs## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=11 +datum=NAD83 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["North American Datum 1983",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6269]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 11N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-117,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16011]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]spplot(cache)
Crop Cache raster to North Logan:
# North Logan Crop:
cache_NL <- crop(cache, NL.sh, mask = T)
cache_NL@extent## class : Extent
## xmin : 928274.6
## xmax : 934574.6
## ymin : 4635899
## ymax : 4640249spplot(cache_NL)
# North Logan Mask
nl_clean <- mask(cache_NL, mask = NL.sh)
spplot(nl_clean)
# Logan City Crop:
cache_L <- crop(cache, L.sh, mask = T)
cache_L@extent## class : Extent
## xmin : 923924.6
## xmax : 934124.6
## ymin : 4626389
## ymax : 4640129spplot(cache_L)
# Logan City Mask:
cache_L_mask <- mask(cache_L, mask = L.sh)
spplot(cache_L_mask)
Q4: There are several categories that indicate “urban” land:
* 121 - Developed/Open Space
* 122 = Developed/Low Intensity
* 123 = Developed/Medium Intensity
* 124 = Developed/High Intensity
Using your new dataset cropped to the extent of the Logan
municipality, reclassify pixels with these values with the value of
999. How many pixels are in this new reclassified
999 urban category?
There are 6,579 pixels reclassified as Urban.
There are 30,450 pixels total in the new urban_cache raster.
# North Logan Reclassify:
urban_cache <- reclassify(nl_clean, c(121, 124, 999))
spplot(urban_cache)
length(urban_cache == 999)## [1] 30450# Logan Reclassify:
logan_urban <- reclassify(cache_L_mask, c(121, 124, 999))
spplot(logan_urban)
length(logan_urban == 999)## [1] 155720Q5: Mask this new urban category, plot the data, and upload this plot on Canvas.
urban <- mask(urban_cache, urban_cache == 999, maskvalue = F)
spplot(urban)
Q6: Resample to 90 meter resolution. How many pixels are in your resampled dataset?
urban_90m <- aggregate(urban_cache, fact =3, fun = modal, na.rm = T)
urban_cache_90m <- resample(urban_cache, urban_90m, method = "ngb")
spplot(urban_cache_90m)
length(urban_cache_90m)## [1] 3430