Module 12: Spatial data
Katya Podkovyroff Lewis
2022-06-22
Spatial data exercises
library(sf)## Linking to GEOS 3.9.1, GDAL 3.4.0, PROJ 8.1.1; sf_use_s2() is TRUElibrary(ncdf4)
library(RColorBrewer)
library(viridis)## Loading required package: viridisLitelibrary(ggplot2)
library(viridisLite)
library(ggthemes)
library(raster)## Loading required package: spExercise 1
In the zipfile countries.zip is a global shapefile with various socio-economic indicators for different countries. Load this file into R and make plots of any two of the following variables (the variable names are given in brackets). Try different color scales and transformations of the data to get the most informative maps
countries <- st_read(dsn = "./countries.shp",
layer = "countries",
drivers = "ESRI Shapefile")## options: ESRI Shapefile
## Reading layer `countries' from data source
## `/Volumes/KATYA4 2022/Intro to R/Module 12/countries.shp' using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## CRS: NAcountries <- st_set_crs(countries, 4326)
st_crs(countries) <- 4326
countries <- read_sf("./countries.shp")
class(countries)## [1] "sf" "tbl_df" "tbl" "data.frame"countries## Simple feature collection with 177 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## CRS: NA
## # A tibble: 177 × 65
## SP_ID scalerank featurecla labelrank sovereignt sov_a3 adm0_dif level type
## <chr> <int> <chr> <dbl> <chr> <chr> <dbl> <dbl> <chr>
## 1 0 1 Admin-0 cou… 3 Afghanist… AFG 0 2 Sove…
## 2 1 1 Admin-0 cou… 3 Angola AGO 0 2 Sove…
## 3 2 1 Admin-0 cou… 6 Albania ALB 0 2 Sove…
## 4 3 1 Admin-0 cou… 4 United Ar… ARE 0 2 Sove…
## 5 4 1 Admin-0 cou… 2 Argentina ARG 0 2 Sove…
## 6 5 1 Admin-0 cou… 6 Armenia ARM 0 2 Sove…
## 7 6 1 Admin-0 cou… 4 Antarctica ATA 0 2 Inde…
## 8 7 3 Admin-0 cou… 6 France FR1 1 2 Depe…
## 9 8 1 Admin-0 cou… 2 Australia AU1 1 2 Coun…
## 10 9 1 Admin-0 cou… 4 Austria AUT 0 2 Sove…
## # … with 167 more rows, and 56 more variables: admin <chr>, adm0_a3 <chr>,
## # geou_dif <dbl>, geounit <chr>, gu_a3 <chr>, su_dif <dbl>, subunit <chr>,
## # su_a3 <chr>, brk_diff <dbl>, name <chr>, name_long <chr>, brk_a3 <chr>,
## # brk_name <chr>, brk_group <chr>, abbrev <chr>, postal <chr>,
## # formal_en <chr>, formal_fr <chr>, note_adm0 <chr>, note_brk <chr>,
## # name_sort <chr>, name_alt <chr>, mapcolor7 <dbl>, mapcolor8 <dbl>,
## # mapcolor9 <dbl>, mapcolor13 <dbl>, pop_est <dbl>, gdp_md_est <dbl>, …countries2 <- countries[ , c("gdp_md_est", "pop_est", "income_grp")]
names(countries2)## [1] "gdp_md_est" "pop_est" "income_grp" "geometry"geometry <- countries2$geometry
Median_GDP <- subset(countries, select = c(gdp_md_est))
Median_GDP <- st_set_crs(Median_GDP, 4326)
st_crs(Median_GDP) <- 4326
names(Median_GDP)## [1] "gdp_md_est" "geometry"plot(Median_GDP, main = "Median GDP as a socio-economic factor")
ggplot() + geom_sf(data = Median_GDP, aes(fill = gdp_md_est)) + theme_bw() + scale_fill_viridis(option = "viridis") + ggtitle("Median GDP as a socio-economic factor") 
Population <- subset(countries, select = c(pop_est))
Population <- st_set_crs(Population, 4326)
st_crs(Population) <- 4326
names(Population)## [1] "pop_est" "geometry"plot(Population, main = "Population as a socio-economic factor")
ggplot() + geom_sf(data = Population, aes(fill = pop_est)) + theme_bw() + scale_fill_viridis(option = "magma") + ggtitle ("Population as a socio-economic factor")
Exercise 2
The NetCDF file air.mon.mean.nc contains air temperature data from the NCAR NCEP project, but as mean monthly temperature for every year between 1948 and the present. Using the functions introduced in this lab, read in the data and make a time series of globally average temperature from January 1948 to present.
air_temp <- raster("./air.mon.mean.nc")
air_temp## class : RasterLayer
## band : 1 (of 819 bands)
## dimensions : 73, 144, 10512 (nrow, ncol, ncell)
## resolution : 2.5, 2.5 (x, y)
## extent : -1.25, 358.75, -91.25, 91.25 (xmin, xmax, ymin, ymax)
## crs : NA
## source : air.mon.mean.nc
## names : Monthly.Mean.Air.Temperature.at.sigma.level.0.995
## z-value : 1948-01-01
## zvar : airhead(air_temp)## 1 2 3 4 5 6
## 1 -34.926773 -34.9267731 -34.9267731 -34.9267731 -34.9267731 -34.9267731
## 2 -35.139351 -35.1296730 -35.1274185 -35.1261215 -35.1306419 -35.1370926
## 3 -34.352573 -34.0422592 -33.7687073 -33.5322571 -33.3509674 -33.2209702
## 4 -27.584831 -26.4677372 -25.5093536 -24.7706451 -24.2719364 -24.0554848
## 5 -16.149351 -14.5541849 -13.3174152 -12.5203142 -12.2361193 -12.5028973
## 6 -8.645802 -7.3235450 -6.5003181 -6.2048311 -6.4119301 -7.1280575
## 7 -4.576451 -4.0216079 -3.9496732 -4.1796761 -4.6077394 -5.1299934
## 8 -2.550969 -2.2787063 -2.0212893 -1.7567741 -1.6116098 -1.6719342
## 9 -1.189029 -0.6406442 0.2729049 0.9612928 0.8025867 -0.4987066
## 10 1.213550 2.1406481 3.0048442 2.5832326 0.2051657 -3.7277381
## 7 8 9 10 11 12 13
## 1 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773
## 2 -35.150318 -35.167740 -35.186443 -35.200966 -35.222256 -35.250969 -35.269035
## 3 -33.147419 -33.132900 -33.175800 -33.263866 -33.388386 -33.543549 -33.716129
## 4 -24.115164 -24.449350 -25.026453 -25.784836 -26.674837 -27.616776 -28.547421
## 5 -13.329993 -14.669668 -16.440956 -18.494513 -20.674839 -22.792583 -24.674839
## 6 -8.316122 -9.928704 -11.846441 -13.944183 -16.048700 -17.958382 -19.501608
## 7 -5.695477 -6.259995 -6.791286 -7.288380 -7.760322 -8.200317 -8.598706
## 8 -1.930642 -2.255158 -2.480001 -2.512256 -2.377096 -2.210000 -2.170645
## 9 -2.665803 -5.057095 -6.994190 -8.064187 -8.249349 -7.773220 -6.968704
## 10 -8.001933 -11.348704 -13.188704 -13.661285 -13.296121 -12.480316 -11.320641
## 14 15 16 17 18 19 20
## 1 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773 -34.926773
## 2 -35.286125 -35.310318 -35.332253 -35.344841 -35.357414 -35.363228 -35.379353
## 3 -33.899029 -34.072254 -34.231934 -34.378067 -34.496452 -34.589352 -34.653217
## 4 -29.400324 -30.119997 -30.684191 -31.063866 -31.276777 -31.334843 -31.286129
## 5 -26.157419 -27.152584 -27.606449 -27.543224 -27.073229 -26.317421 -25.457417
## 6 -20.503550 -20.886127 -20.637737 -19.868063 -18.755480 -17.562899 -16.574831
## 7 -8.910962 -9.068057 -9.039024 -8.835799 -8.591928 -8.511604 -8.870317
## 8 -2.346128 -2.696769 -3.115805 -3.501938 -3.882899 -4.434191 -5.439990
## 9 -6.077092 -5.263224 -4.589353 -4.133225 -3.971932 -4.224837 -5.044835
## 10 -9.846122 -8.238381 -6.939027 -6.408057 -6.839351 -8.017739 -9.500640air_temp = rotate(air_temp)
my.pal <- brewer.pal(n = 9, name = "PiYG")
plot(air_temp,
main = "Average Global Temperatures January 1948-Present",
col = my.pal)
nc_countries <- st_read("./ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp",
quiet = TRUE)
plot(st_geometry(nc_countries), add = TRUE)