In Robinlovelace/geocompr: Geocomputation with R

library(sf)
library(terra)
library(spData)

E1. Create a new object called nz_wgs by transforming nz object into the WGS84 CRS.

Create an object of class crs for both and use this to query their CRSs.
With reference to the bounding box of each object, what units does each CRS use?
Remove the CRS from nz_wgs and plot the result: what is wrong with this map of New Zealand and why?

st_crs(nz)
nz_wgs = st_transform(nz, "EPSG:4326")
nz_crs = st_crs(nz)
nz_wgs_crs = st_crs(nz_wgs)
nz_crs$epsg
nz_wgs_crs$epsg
st_bbox(nz)
st_bbox(nz_wgs)
nz_wgs_NULL_crs = st_set_crs(nz_wgs, NA)
nz_27700 = st_transform(nz_wgs, "EPSG:27700")
par(mfrow = c(1, 3))
plot(st_geometry(nz))
plot(st_geometry(nz_wgs))
plot(st_geometry(nz_wgs_NULL_crs))
# answer: it is fatter in the East-West direction
# because New Zealand is close to the South Pole and meridians converge there
plot(st_geometry(nz_27700))
par(mfrow = c(1, 1))

E2. Transform the world dataset to the transverse Mercator projection ("+proj=tmerc") and plot the result. What has changed and why? Try to transform it back into WGS 84 and plot the new object. Why does the new object differ from the original one?

# see https://github.com/r-spatial/sf/issues/509
world_tmerc = st_transform(world, "+proj=tmerc")
plot(st_geometry(world_tmerc))
world_4326 = st_transform(world_tmerc, "EPSG:4326")
plot(st_geometry(world_4326))

E3. Transform the continuous raster (con_raster) into NAD83 / UTM zone 12N using the nearest neighbor interpolation method. What has changed? How does it influence the results?

con_raster = rast(system.file("raster/srtm.tif", package = "spDataLarge"))
con_raster_utm12n = project(con_raster, "EPSG:32612", method = "near")
con_raster_utm12n

plot(con_raster)
plot(con_raster_utm12n)

E4. Transform the categorical raster (cat_raster) into WGS 84 using the bilinear interpolation method. What has changed? How does it influence the results?

cat_raster = rast(system.file("raster/nlcd.tif", package = "spDataLarge"))
cat_raster_wgs84 = project(cat_raster, "EPSG:4326", method = "bilinear")
cat_raster_wgs84

plot(cat_raster)
plot(cat_raster_wgs84)