Chapter 6
Harold Nelson
5/30/2021
Packages
# install.packages("rasterVis")
library(sf) # Classes a methods for spatial data - simple features
library(raster) # Classes and methods for raster data
library(rgdal) # functions for spatial data input/output
library(RColorBrewer) # Creates nice color schemes
library(rasterVis) # Advanced plotting functions for raster objects
library(ggplot2)
library(tmap)
library(tmaptools)
Get lst data
lst <- raster("MOD11A2_2017-07-12.LST_Day_1km.tif")
class(lst)
## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"
## class : RasterLayer
## dimensions : 1110, 3902, 4331220 (nrow, ncol, ncell)
## resolution : 0.009009009, 0.009009009 (x, y)
## extent : -104.4326, -69.27943, 30, 40 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : MOD11A2_2017-07-12.LST_Day_1km.tif
## names : MOD11A2_2017.07.12.LST_Day_1km
## values : 0, 65535 (min, max)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (2.31% of all cells)
## MOD11A2_2017.07.12.LST_Day_1km
## Min. 0
## 1st Qu. 0
## Median 14981
## 3rd Qu. 15148
## Max. 16050
## NA's 0
Convert
lst[lst == 0] <- NA
lst_c <- lst * 0.02 - 273.16
hist(lst_c)
Get US States
Then use tmap to include a raster.
ct = st_read("conterminous_us_states.shp")
## Reading layer `conterminous_us_states' from data source
## `/Users/haroldnelson/Dropbox/Wimberly/gdswr_chapter6_data/conterminous_us_states.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 49 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.7328 ymin: 24.95638 xmax: -66.96927 ymax: 49.37173
## CRS: NA
tm_shape(ct) +
tm_polygons() +
tm_shape(lst_c) +
tm_raster()
## Warning: Currect projection of shape ct unknown. Long-lat (WGS84) is assumed.
## stars object downsampled to 1875 by 533 cells. See tm_shape manual (argument raster.downsample)
Apparently this is a strange wedge-shaped region.
Get a Brick
## NLDAS_FORA0125_M.A201207.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201307.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201407.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201507.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201607.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201707.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201207.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201307.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201407.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201507.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201607.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
## NLDAS_FORA0125_M.A201707.002.grb has GDAL driver GRIB
## and has 224 rows and 464 columns
Scatterplot
Do a scatterplot of two layers.
plot(temp2012, temp2013, xlab="2012", ylab = "2013")
## Warning in .local(x, y, ...): plot used a sample of 95.9% of the cells. You can
## use "maxpixels" to increase the sample)
Make a Brick
tempbrick <- brick(temp2012, temp2013, temp2014,
temp2015, temp2016, temp2017)
names(tempbrick) <- c("July.2012", "July.2013", "July.2014",
"July.2015", "July.2016", "July.2017")
tempbrick
## class : RasterBrick
## dimensions : 224, 464, 103936, 6 (nrow, ncol, ncell, nlayers)
## resolution : 0.125, 0.125 (x, y)
## extent : -125.0005, -67.0005, 25.0005, 53.0005 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +R=6371200 +no_defs
## source : memory
## names : July.2012, July.2013, July.2014, July.2015, July.2016, July.2017
## min values : 4.150000, 3.180000, 3.320000, 4.010000, 3.090000, 5.270039
## max values : 37.35000, 38.95000, 39.17000, 36.63000, 38.83000, 40.18004
rasterdf
This is not needed if you use tmap.
rasterdf <- function(x, aggregate = 1) {
resampleFactor <- aggregate
inputRaster <- x
inCols <- ncol(inputRaster)
inRows <- nrow(inputRaster)
# Compute numbers of columns and rows in the new raster for mapping
resampledRaster <- raster(ncol=(inCols / resampleFactor),
nrow=(inRows / resampleFactor))
# Match to the extent of the original raster
extent(resampledRaster) <- extent(inputRaster)
# Resample data on the new raster
y <- resample(inputRaster,resampledRaster,method='ngb')
# Extract cell coordinates into a data frame
coords <- xyFromCell(y, seq_len(ncell(y)))
# Extract layer names
dat <- stack(as.data.frame(getValues(y)))
# Add names - 'value' for data, 'variable' to indicate different raster layers
# in a stack
names(dat) <- c('value', 'variable')
dat <- cbind(coords, dat)
dat
}
tempbrick_df <- rasterdf(tempbrick)
str(tempbrick)
## Formal class 'RasterBrick' [package "raster"] with 12 slots
## ..@ file :Formal class '.RasterFile' [package "raster"] with 13 slots
## .. .. ..@ name : chr ""
## .. .. ..@ datanotation: chr "FLT4S"
## .. .. ..@ byteorder : chr "little"
## .. .. ..@ nodatavalue : num -Inf
## .. .. ..@ NAchanged : logi FALSE
## .. .. ..@ nbands : int 1
## .. .. ..@ bandorder : chr "BIL"
## .. .. ..@ offset : int 0
## .. .. ..@ toptobottom : logi TRUE
## .. .. ..@ blockrows : int 0
## .. .. ..@ blockcols : int 0
## .. .. ..@ driver : chr ""
## .. .. ..@ open : logi FALSE
## ..@ data :Formal class '.MultipleRasterData' [package "raster"] with 14 slots
## .. .. ..@ values : num [1:103936, 1:6] 14.6 14.1 13.3 12.8 13.5 ...
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : NULL
## .. .. .. .. ..$ : chr [1:6] "band1.1" "band1.2" "band1.3" "band1.4" ...
## .. .. ..@ offset : num 0
## .. .. ..@ gain : num 1
## .. .. ..@ inmemory : logi TRUE
## .. .. ..@ fromdisk : logi FALSE
## .. .. ..@ nlayers : int 6
## .. .. ..@ dropped : NULL
## .. .. ..@ isfactor : logi FALSE
## .. .. ..@ attributes: list()
## .. .. ..@ haveminmax: logi TRUE
## .. .. ..@ min : num [1:6] 4.15 3.18 3.32 4.01 3.09 ...
## .. .. ..@ max : num [1:6] 37.4 39 39.2 36.6 38.8 ...
## .. .. ..@ unit : chr ""
## .. .. ..@ names : chr [1:6] "July.2012" "July.2013" "July.2014" "July.2015" ...
## ..@ legend :Formal class '.RasterLegend' [package "raster"] with 5 slots
## .. .. ..@ type : chr(0)
## .. .. ..@ values : logi(0)
## .. .. ..@ color : logi(0)
## .. .. ..@ names : logi(0)
## .. .. ..@ colortable: logi(0)
## ..@ title : chr(0)
## ..@ extent :Formal class 'Extent' [package "raster"] with 4 slots
## .. .. ..@ xmin: num -125
## .. .. ..@ xmax: num -67
## .. .. ..@ ymin: num 25
## .. .. ..@ ymax: num 53
## ..@ rotated : logi FALSE
## ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
## .. .. ..@ geotrans: num(0)
## .. .. ..@ transfun:function ()
## ..@ ncols : int 464
## ..@ nrows : int 224
## ..@ crs :Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +R=6371200 +no_defs"
## .. .. ..$ comment: chr "GEOGCRS[\"Coordinate System imported from GRIB file\",\n DATUM[\"unnamed\",\n ELLIPSOID[\"Sphere\",63"| __truncated__
## ..@ history : list()
## ..@ z : list()
tmap
The preceding step is easier with tmap.
st_crs(ct) = 4326
tempbrick = mask(tempbrick,ct)
tm_shape(tempbrick) +
tm_raster() +
tm_shape(ct) +
tm_borders()
Mean temperature
Get the mean July temperature for all six years and map with tmap.
meantemp <- mean(tempbrick)
tm_shape(meantemp) +
tm_raster() +
tm_shape(ct) +
tm_borders() +
tm_layout(legend.outside = TRUE,
legend.outside.position = c("left","bottom"),
title = "Celsius",
main.title = "Mean July Temp",)
Anomalies
tempanom <- tempbrick - meantemp
names(tempanom) <- names(tempbrick)
tm_shape(tempanom) +
tm_raster() +
tm_shape(ct) +
tm_borders()
## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
