Extents & discritizations
Mike Johnson
Lynker, NOAA-AffiliateSource:
vignettes/region.Rmd
region.Rmd
library(AOI)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, unionBasic Usage
The core functions described here are aoi_ext which
builds a specific extent around a given POINT /
coordinate set and discritize which
discritizes an extent/object.
To start, lets build a 1000 meter extent around the location called “Fort Collins”:
aoi_ext(geo = "Fort Collins", wh = 1000)
#> xmin ymin xmax ymax
#> -105.08967 40.57283 -105.06605 40.59084From this, we can see a extent was created!
Add bbox
As with AOI::geocode, the returned object can be
represented as a POLYGON by adding
bbox = TRUE. This is helpful for mapping:
Modify both diminisons
Providing a single value to the wh argument generates a
square extent. To find a rectangular extent, a two value vector can be
passed:
Use XY
We have relied on providing a place that requires
AOI::geocode to be run. Instead, we can provide a direct
coordinate set as an xy vector:
Units
By default all measures are made in the default AOI unit:
(AOI:::default_units)
#> [1] "meter"We can override these by defining the units argument in
the aoi_ext signiture. Here we see that a 1km and 1000m wh argument
produce identical results:
Area discritizations
In many applications, getting the extent is only part of the
challenge. The other is defining a discretizion of a region. This can be
done with the AOI::discritize function:
aoi_ext("Fort Collins", wh = c(10000, 20000)) %>%
discritize()
#> $extent
#> xmin ymin xmax ymax
#> -105.19630 40.40167 -104.95943 40.76188
#>
#> $dimension
#> [1] 337 512
#>
#> $projection
#> [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"Diminsions
By default, the function uses the AOI default discritization:
(AOI:::default_dim)
#> [1] 512This can be changed using the dim argument:
aoi_ext("Fort Collins", wh = c(10000, 20000)) %>%
discritize(dim = 1024)
#> $extent
#> xmin ymin xmax ymax
#> -105.19630 40.40167 -104.95943 40.76188
#>
#> $dimension
#> [1] 673 1024
#>
#> $projection
#> [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"CRS
The default expected input CRS is EPSG:4326. This is also the default output. Both can be changed depending on the desired input and output.
- If a CRS can be extracted from the input, then it is used.
- If a CRS cannot be extracted from the input, then EPSG:4326 is assumed
- The assumed CRS can be overridden with the
in_crsargument
# CRS inferred from input and changed to 5070
aoi_ext("Fort Collins", wh = c(20000, 10000)) %>%
discritize(out_crs = 5070)
#> $extent
#> xmin ymin xmax ymax
#> -781559.4 1976585.3 -740107.0 2000454.5
#>
#> $dimension
#> [1] 512 889
#>
#> $projection
#> [1] "PROJCRS[\"NAD83 / Conus Albers\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"Conus Albers\",\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 (X)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (Y)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Data analysis and small scale data presentation for contiguous lower 48 states.\"],\n AREA[\"United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming.\"],\n BBOX[24.41,-124.79,49.38,-66.91]],\n ID[\"EPSG\",5070]]"
# CRS inferred
t1 = aoi_ext("Fort Collins", wh = c(20000, 10000), crs = 5070, bbox = TRUE) %>%
discritize(out_crs = 5070)
t1$extent
#> xmin ymin xmax ymax
#> -781506.2 1976590.7 -740107.0 2000454.0
# Bad, assumed EPSG
t2 = c(-781506.2, 1976590.7, -740107.0, 2000454.0) %>%
discritize(out_crs = 5070)
#> Warning: Assuming the input is in CRS 4326
t2$extent
#> xmin ymin xmax ymax
#> NA NA NA NA
# Explicit CRS named
t3 = c(-781506.2, 1976590.7, -740107.0, 2000454.0) %>%
discritize(in_crs = 5070, out_crs = 5070)
t3$extent
#> xmin ymin xmax ymax
#> -781506.2 -740107.0 1976590.7 2000454.0SpatRaster
If you want a gridded object for future processing, you can set
spatrast = TRUE.
aoi_ext("Fort Collins", wh = c(20000, 10000)) %>%
discritize(spatrast = TRUE)
#> class : SpatRaster
#> dimensions : 1345, 512, 1 (nrow, ncol, nlyr)
#> resolution : 0.0009239945, 0.0001339108 (x, y)
#> extent : -105.3144, -104.8413, 40.49154, 40.67165 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)Passing Other objects.
Defining and extent for a POINT and wh
argument is only one option that can be discritized. Others include all
AOI outputs, sf objects and both SpatVect and SpatRast objects.
AOI BBOX returns
geocode(c("Fort Collins", "Boulder"), bbox = TRUE) %>%
discritize(spatrast = TRUE, dim = c(10,10))
#> class : SpatRaster
#> dimensions : 10, 4, 1 (nrow, ncol, nlyr)
#> resolution : 0.050655, 0.056573 (x, y)
#> extent : -105.2792, -105.0766, 40.01574, 40.58147 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)AOI fiat boundaries
(really this applies to any sf object!)
aoi_get(state = "OR") %>%
discritize(spatrast = TRUE, dim = c(5,5))
#> class : SpatRaster
#> dimensions : 10, 5, 1 (nrow, ncol, nlyr)
#> resolution : 1.648007, 0.4307017 (x, y)
#> extent : -124.7035, -116.4635, 41.99208, 46.2991 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)Other Spat* object
(r <- terra::rast(system.file("ex/elev.tif", package="terra")) )
#> class : SpatRaster
#> dimensions : 90, 95, 1 (nrow, ncol, nlyr)
#> resolution : 0.008333333, 0.008333333 (x, y)
#> extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : elev.tif
#> name : elevation
#> min value : 141
#> max value : 547
# SpatRast
discritize(r, spatrast = TRUE)
#> class : SpatRaster
#> dimensions : 540, 512, 1 (nrow, ncol, nlyr)
#> resolution : 0.001546224, 0.001388889 (x, y)
#> extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
# SpatVect
terra::vect(system.file("ex/lux.shp", package="terra")) %>%
discritize(spatrast = TRUE)
#> class : SpatRaster
#> dimensions : 547, 512, 1 (nrow, ncol, nlyr)
#> resolution : 0.001531469, 0.001341525 (x, y)
#> extent : 5.74414, 6.528252, 49.44781, 50.18162 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)