Neon canopy height from space (solution)

Exercise

Output solution


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE

terra 1.8.70


Attaching package: 'tidyterra'

The following object is masked from 'package:stats':

    filter

1 Map the digital terrain model for SJER using the viridis color ramp.

2 Create and map the canopy height model for SJER using raster math (chm = dsm - dtm) and the viridis color ramp.

3 Create a map that shows the SJER boundary and the plot locations colored by plot type.

Transform the plot data to have the same CRS as the CHM. Show CRS for both objects. Create a map that shows the canopy height model from (3) with the plot locations on top.

Coordinate Reference System:
  User input: WGS 84 / UTM zone 11N 
  wkt:
PROJCRS["WGS 84 / UTM zone 11N",
    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["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]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 120°W and 114°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; British Columbia (BC); Northwest Territories (NWT); Nunavut. Mexico. United States (USA)."],
        BBOX[0,-120,84,-114]],
    ID["EPSG",32611]]

Coordinate Reference System:
  User input: WGS 84 / UTM zone 11N 
  wkt:
PROJCRS["WGS 84 / UTM zone 11N",
    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["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]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 120°W and 114°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Alberta; British Columbia (BC); Northwest Territories (NWT); Nunavut. Mexico. United States (USA)."],
        BBOX[0,-120,84,-114]],
    ID["EPSG",32611]]

Extract the canopy heights at each plot location for SJER and display the values.

[1] 18.913757 23.948151  1.986877  2.183136 28.985016  3.506866  2.201233

Building on (5) create a version of your code that extracts the canopy heights and includes them with the data in plots_sjer_utm. You can do this using either mutate (to add the results of (5) to plots_sjer_utm or bind put the data together in extract and st_as_sf to convert the resulting terra object into a simple features object. Display the full simple features object. Make sure that the column name for the canopy heights is canopy_height. The detailed display will vary depending on your approach.

Simple feature collection with 7 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 256238.5 ymin: 4110088 xmax: 257465.5 ymax: 4111372
Projected CRS: WGS 84 / UTM zone 11N
  plot_id   plot_type canopy_heights                 geometry
1       1       Tower      18.913757 POINT (257465.5 4111372)
2       2       Tower      23.948151 POINT (256238.5 4110270)
3       3       Tower       1.986877 POINT (256798.3 4110296)
4       4 Distributed       2.183136 POINT (256737.7 4110949)
5       5 Distributed      28.985016 POINT (257358.5 4110450)
6       6 Distributed       3.506866 POINT (256254.5 4110088)
7       7 Distributed       2.201233 POINT (256754.5 4110274)

Create a map that shows the SJER boundary and the plot locations colored by the canopy height values.

Create a map that shows the canopy height model raster, but in cm rather than m (i.e., multiply the canopy height model by 100).

Create a map that shows the digital terrain model (DTM) raster, the plot locations, and the SJER boundary, using transparency as needed to allow all three layers to be seen. Remember all three layers will need to have the same CRS.

Conduct an analysis of the relationship between elevation and canopy height at the SJER plots. Using a 50m buffer, extract the mean elevations (i.e., the values from the digital terrain model) and the canopy heights at each plot location for SJER and add to the spatial plots data to produce a simple features object that includes both the average elevations (in a 50 m buffer) and the canopy heights (in a 50 m buffer). Then make a scatter plot showing the relationship between elevation and canopy height using this data. Color the points by plot type and fit a single smooth curve through all of the points. Finally, use dplyr to calculate the average canopy height and average elevation for the two different plot types.

`geom_smooth()` using formula = 'y ~ x'

Simple feature collection with 2 features and 3 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: 256238.5 ymin: 4110088 xmax: 257465.5 ymax: 4111372
Projected CRS: WGS 84 / UTM zone 11N
# A tibble: 2 × 4
  plot_type   elevation canopy_height                                   geometry
  <chr>           <dbl>         <dbl>                           <MULTIPOINT [m]>
1 Distributed      383.          2.06 ((256254.5 4110088), (256737.7 4110949), …
2 Tower            386.          2.63 ((256238.5 4110270), (256798.3 4110296), …