1. Load each dataset
2. Check geometry validity
3. Align CRS
4. Run Correlation
5. Print ResultsAssignment 3 Solutions: Coordinates and Geometries
1. Write out the pseudocode that you would use to set up an analysis of the spatial correlations between chronic asthma risk, exposure to PM2.5, and wildfire. You don’t have to write functions or any actual code. Just write the steps and insert named code blocks for each step.
This one is probably a little tricky if you haven’t taken the time to check out the attributes of the data (which you should always do). That said, some pretty generic steps would be:
There are two key steps here, that you’ll repeat for any/all spatial analyses that you do: 1) checking for valid geometries and 2) making sure the data are aligned in a sensible CRS. I can add a code chunk for each now.
1. Load each dataset2. Check geometry validity3. Align CRS4. Run Correlation5. Print Results2. Read in the cdc_nw.shp, pm_nw.shp, and wildfire_hazard_agg.tif files and print the coordinate reference system for each object. Do they match?
Here I’m going to combine the
loadportion of my pseudocode with thevaliditysince I can do that without creating additional object. I use thestr()function to get a sense for what the data looks like and to understand what data classes I’m working with. Then, I use theall()function to make sure that all of the results ofst_is_valid()are true. I don’t need to do that with the raster file as the geometry is implicit which means that it has to be topologically valid (this doesn’t mean that the numbers are accurate, it just means that the dataset conforms to the data modelRexpects). Then I’ll add another code to check theCRSof the different objects.
library(sf)
library(terra)
library(tidyverse)
cdc.nw <- read_sf("data/opt/data/2023/assignment03/cdc_nw.shp")
str(cdc.nw)
all(st_is_valid(cdc.nw))
pm.nw <- read_sf("data/opt/data/2023/assignment03/pm_nw.shp")
str(pm.nw)
all(st_is_valid(pm.nw))
wildfire.haz <- rast("data/opt/data/2023/assignment03/wildfire_hazard_agg.tif")
str(wildfire.haz)Now that I’ve gotten the data into my environment, I need to make sure that the CRS are aligned. I’ll demonstrate that with a few different approaches. You can use the logical
==or theidenticalfunction to check, but remember that these functions are not specific to spatial objects, they evaluate things very literally. So even if the CRS is the same, ifst_crsreturns the CRS in one format (WKT) andcrsreturns it in another, you’ll getFALSEeven if they are actually the same CRS - pay attention to that. You’ll notice that they aren’t identical; we’ll deal with that in the next question.
st_crs(cdc.nw)Coordinate Reference System:
User input: NAD83
wkt:
GEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]]st_crs(pm.nw)Coordinate Reference System:
User input: WGS 84 / Pseudo-Mercator
wkt:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]crs(wildfire.haz)[1] "PROJCRS[\"unnamed\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"Albers Equal Area\",\n METHOD[\"Albers Equal Area\",\n ID[\"EPSG\",9822]],\n PARAMETER[\"Latitude of false origin\",23,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",-96,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",29.5,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",45.5,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"identical(st_crs(cdc.nw), st_crs(pm.nw))[1] FALSEst_crs(cdc.nw) == st_crs(pm.nw)[1] FALSE3. Re-project the cdc_nw.shp and pm_nw.shp shapefiles so that they have the same CRS as the wildfire_hazard_agg.tfi file. Verify that all the files have the same projection.
Now we’ll use
st_transformto get the two shapefiles aligned with the raster (because we generally want to avoid projecting rasters if we can). We can then use the same steps above to see if they’re aligned. Note that I’m using theterra::crs()function to make sure that the output is printed in exactly the same format
cdc.nw.proj <- cdc.nw %>% st_transform(., crs=crs(wildfire.haz))
pm.nw.proj <- pm.nw %>% st_transform(., crs=crs(wildfire.haz))
identical(crs(cdc.nw.proj), crs(wildfire.haz))[1] TRUEidentical(crs(pm.nw.proj), crs(wildfire.haz))[1] TRUE4. How does reprojecting change the coordinates of the bounding box for the two shapefiles? Show your code
Now we just want to look at the bounding box of the data before and after it was projected. We can do this using
st_bbox. One of the most obvious changes is that the units forcdc.nwhave changed from degrees to meters (as evidenced by the much larger numbers). For thepm.nwobject we can see that the raw coordinates indicate a shift to the west; however, because the origin for this crs has also changed, the states still show up in the correct place.
st_bbox(cdc.nw) xmin ymin xmax ymax
-124.74918 41.98818 -111.04349 49.00232 st_bbox(cdc.nw.proj) xmin ymin xmax ymax
-2295337 2208890 -1189292 3177425 st_bbox(pm.nw) xmin ymin xmax ymax
-13898126 5159210 -12361307 6275276 st_bbox(pm.nw.proj) xmin ymin xmax ymax
-2300791 2208891 -1189293 3177426 5. What class of geometry does the pm_nw.shp have (show your code)? Now filter the pm_nw.shp file so that only the records from Ada County, Idaho are showing. Find the record with the lowest value for PM25. How many coordinates are associated with that geometry?
This one was probably a little tricky. First, to check the geometry type, we use
st_geometry_typesettingby_geometrytoFALSEmeans we get the geometry type for the entire object instead of each observation. We then use a series offiltercommands to get the records from Idaho and Ada county. Once we’ve narrowed the data to our correct region, we canfilteragain to find the row with the minimum value of PM25 (note that we have to setna.rm=TRUEso that we ignore the NA values). Then we just take the number of rows (nrow) of the result ofst_coordinatesto get the number of coordinates associated with that geometry.
st_geometry_type(pm.nw, by_geometry = FALSE)[1] MULTIPOLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLEada.pm <- pm.nw %>%
filter(STATE_NAME=="Idaho" & CNTY_NAME=="Ada") %>%
filter(PM25 == min(PM25, na.rm = TRUE))
ada.pmSimple feature collection with 1 feature and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -12935300 ymin: 5329192 xmax: -12910260 ymax: 5402433
Projected CRS: WGS 84 / Pseudo-Mercator
# A tibble: 1 × 4
STATE_NAME CNTY_NAME PM25 geometry
* <chr> <chr> <dbl> <MULTIPOLYGON [m]>
1 Idaho Ada 6.68 (((-12935301 5391002, -12934885 5391290, -12934526…nrow(st_coordinates(ada.pm))[1] 1394