Spatial Data and Cartography (Part 2)

Lecture 18

Dr. Colin Rundel

Plotting

Example Data - NC SIDS

( nc = read_sf(system.file("shape/nc.shp", package="sf"), quiet = TRUE) |> 
    select(-(AREA:CNTY_ID), -(FIPS:CRESS_ID)))

Simple feature collection with 100 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
# A tibble: 100 × 8
   NAME        BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
   <chr>       <dbl> <dbl>   <dbl> <dbl> <dbl>   <dbl>
 1 Ashe         1091     1      10  1364     0      19
 2 Alleghany     487     0      10   542     3      12
 3 Surry        3188     5     208  3616     6     260
 4 Currituck     508     1     123   830     2     145
 5 Northampton  1421     9    1066  1606     3    1197
 6 Hertford     1452     7     954  1838     5    1237
 7 Camden        286     0     115   350     2     139
 8 Gates         420     0     254   594     2     371
 9 Warren        968     4     748  1190     2     844
10 Stokes       1612     1     160  2038     5     176
# … with 90 more rows, and 1 more variable:
#   geometry <MULTIPOLYGON [°]>

Base Plots

plot(nc)

Geometry Plot

plot(st_geometry(nc), axes=TRUE)

Graticules

plot(nc[,"SID79"], axes=TRUE)

Graticules

plot(nc[,"SID79"], graticule=st_crs(nc), axes=TRUE)

Graticules (EPSG:3631)

plot(st_transform(nc[,"SID79"], 3631), axes=TRUE)

Graticules (EPSG:3631)

plot(st_transform(nc[,"SID79"], 3631), graticule=st_crs(nc), axes=TRUE)

Graticules (EPSG:3631)

plot(st_transform(nc[,"SID79"], 3631), graticule=st_crs(3631), axes=TRUE)

ggplot2

ggplot(nc) + 
  geom_sf(aes(fill=SID79))

ggplot2 + projections

ggplot(st_transform(nc, 3631)) + 
  geom_sf(aes(fill=SID79 / BIR79))

ggplot2 + viridis

ggplot(st_transform(nc, 3631)) + 
  geom_sf(aes(fill=SID79 / BIR79)) +
  scale_fill_viridis_c()

Example Data - Meuse

data(meuse, meuse.riv, package="sp")
(meuse = st_as_sf(meuse, coords=c("x", "y"), crs=28992) |>
  as_tibble() |> st_as_sf())

Simple feature collection with 155 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 178605 ymin: 329714 xmax: 181390 ymax: 333611
Projected CRS: Amersfoort / RD New
# A tibble: 155 × 13
   cadmium copper  lead  zinc  elev    dist    om ffreq
     <dbl>  <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl> <fct>
 1    11.7     85   299  1022  7.91 0.00136  13.6 1    
 2     8.6     81   277  1141  6.98 0.0122   14   1    
 3     6.5     68   199   640  7.8  0.103    13   1    
 4     2.6     81   116   257  7.66 0.190     8   1    
 5     2.8     48   117   269  7.48 0.277     8.7 1    
 6     3       61   137   281  7.79 0.364     7.8 1    
 7     3.2     31   132   346  8.22 0.190     9.2 1    
 8     2.8     29   150   406  8.49 0.0922    9.5 1    
 9     2.4     37   133   347  8.67 0.185    10.6 1    
10     1.6     24    80   183  9.05 0.310     6.3 1    
# … with 145 more rows, and 5 more variables: soil <fct>,
#   lime <fct>, landuse <fct>, dist.m <dbl>,
#   geometry <POINT [m]>

( meuse_riv = st_polygon(list(meuse.riv)) |>
    st_sfc() |>
    st_set_crs(28992) |>
    st_as_sf()
)

Simple feature collection with 1 feature and 0 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 178304 ymin: 325698.5 xmax: 182331.5 ymax: 337684.8
Projected CRS: Amersfoort / RD New
                               x
1 POLYGON ((182003.7 337678.6...

Meuse

plot(meuse, pch=16)

Layering plots

plot(meuse_riv, col=adjustcolor("lightblue", alpha.f=1), border = NA, 
     axes=TRUE, graticule=st_crs(4326), 
     xlim=st_bbox(meuse)[c(1,3)], ylim=st_bbox(meuse)[c(2,4)])
plot(meuse[,"lead"], pch=16, add=TRUE)

ggplot2

ggplot() +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1)

ggplot2 - axis limits

ggplot() +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(50.95, 50.99)

ggplot2 - axis limits

ggplot() +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(329714, 333611)

ggplot2 - bounding box

ggplot() +
  geom_sf(data=st_sf(meuse_riv), fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(st_bbox(meuse)["ymin"], st_bbox(meuse)["ymax"])

Geometry Predicates

DE-9IM

Spatial predicates

st_within(a,b):

\[ \begin{bmatrix} T & * & F \\ * & * & F \\ * & * & * \end{bmatrix} \]

st_touches(a,b):

\[ \begin{bmatrix} F & T & * \\ * & * & * \\ * & * & * \end{bmatrix} \cup \begin{bmatrix} F & * & * \\ T & * & * \\ * & * & * \end{bmatrix} \cup \begin{bmatrix} F & * & * \\ * & T & * \\ * & * & * \end{bmatrix} \]

Sparse vs Full Results

st_intersects(ncc[20:30,], air) %>% str()

List of 11
 $ : int(0) 
 $ : int(0) 
 $ : int(0) 
 $ : int(0) 
 $ : int(0) 
 $ : int 268
 $ : int 717
 $ : int(0) 
 $ : int(0) 
 $ : int(0) 
 $ : int(0) 
 - attr(*, "predicate")= chr "intersects"
 - attr(*, "region.id")= chr [1:11] "1" "2" "3" "4" ...
 - attr(*, "remove_self")= logi FALSE
 - attr(*, "retain_unique")= logi FALSE
 - attr(*, "ncol")= int 940
 - attr(*, "class")= chr [1:2] "sgbp" "list"

st_intersects(ncc, air, sparse=FALSE) %>% str()

 logi [1:100, 1:940] FALSE FALSE FALSE FALSE FALSE FALSE ...

Examples

• Which counties have an airport?

• Which counties are adjacent to Durham County?

• Which counties have more than 4 neighbors?

ncc  = read_sf("data/gis/nc_counties/", quiet=TRUE)
air = read_sf("data/gis/airports/", quiet=TRUE) |> st_transform(st_crs(ncc))
hwy = read_sf("data/gis/us_interstates/", quiet=TRUE) |> st_transform(st_crs(ncc))

Which counties have an airport?

ncc |>
  mutate(
    airports = st_intersects(ncc, air) |> strip_attrs(),
    n_airports = purrr::map_int(airports, length)
  ) |>
  filter(n_airports > 0)

Simple feature collection with 16 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -82.88783 ymin: 33.92835 xmax: -75.46003 ymax: 36.26165
Geodetic CRS:  NAD83
# A tibble: 16 × 11
     AREA PERIM…¹ COUNT…² STATE COUNTY FIPS  STATE…³ SQUAR…⁴
 *  <dbl>   <dbl>   <dbl> <chr> <chr>  <chr> <chr>     <dbl>
 1 0.107     1.54    2065 NC    Forsy… 37067 37         413.
 2 0.170     1.69    2069 NC    Guilf… 37081 37         658.
 3 0.0996    6.15    2077 NC    Dare … 37055 37         386.
 4 0.221     2.14    2106 NC    Wake … 37183 37         858.
 5 0.169     2.12    2151 NC    Pitt … 37147 37         654.
 6 0.106     1.56    2152 NC    Cataw… 37035 37         413.
 7 0.170     2.37    2153 NC    Bunco… 37021 37         660.
 8 0.143     1.71    2197 NC    Wayne… 37191 37         557.
 9 0.140     1.92    2209 NC    Meckl… 37119 37         547.
10 0.181     2.16    2210 NC    Moore… 37125 37         706.
11 0.0937    1.36    2212 NC    Cabar… 37025 37         365.
12 0.103     1.63    2231 NC    Lenoi… 37107 37         403.
13 0.185     4.15    2236 NC    Crave… 37049 37         722.
14 0.169     2.05    2260 NC    Cumbe… 37051 37         659.
15 0.195     4.53    2314 NC    Onslo… 37133 37         766.
16 0.0519    2.23    2408 NC    New H… 37129 37         205.
# … with 3 more variables: geometry <MULTIPOLYGON [°]>,
#   airports <list>, n_airports <int>, and abbreviated
#   variable names ¹​PERIMETER, ²​COUNTYP010, ³​STATE_FIPS,
#   ⁴​SQUARE_MIL

Which counties are adjacent to Durham County?

ncc |>
  mutate(
    touch_durham = st_touches(ncc, ncc |> filter(COUNTY == "Durham County")) |> strip_attrs(),
    n_touches = map_int(touch_durham, length)
  ) |>
  filter(n_touches > 0)

Simple feature collection with 5 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.55573 ymin: 35.51505 xmax: -78.25503 ymax: 36.54245
Geodetic CRS:  NAD83
# A tibble: 5 × 11
   AREA PERIMETER COUNT…¹ STATE COUNTY FIPS  STATE…² SQUAR…³
* <dbl>     <dbl>   <dbl> <chr> <chr>  <chr> <chr>     <dbl>
1 0.105      1.30    2006 NC    Perso… 37145 37         404.
2 0.139      1.69    2008 NC    Granv… 37077 37         536.
3 0.104      1.30    2074 NC    Orang… 37135 37         401.
4 0.221      2.14    2106 NC    Wake … 37183 37         858.
5 0.183      2.26    2141 NC    Chath… 37037 37         710.
# … with 3 more variables: geometry <MULTIPOLYGON [°]>,
#   touch_durham <list>, n_touches <int>, and abbreviated
#   variable names ¹​COUNTYP010, ²​STATE_FIPS, ³​SQUARE_MIL

Which counties have more than 4 neighbors?

ncc |>
  mutate(
    neighbors = st_touches(ncc) |> strip_attrs(),
    n_neighbors = map_int(neighbors, length)
  ) |>
  ggplot(aes(fill = n_neighbors > 4)) +
    geom_sf()

Geometry Manipulation

Casting

(nc_pts = st_cast(nc, "MULTIPOINT"))

Simple feature collection with 100 features and 7 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
# A tibble: 100 × 8
   NAME        BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
   <chr>       <dbl> <dbl>   <dbl> <dbl> <dbl>   <dbl>
 1 Ashe         1091     1      10  1364     0      19
 2 Alleghany     487     0      10   542     3      12
 3 Surry        3188     5     208  3616     6     260
 4 Currituck     508     1     123   830     2     145
 5 Northampton  1421     9    1066  1606     3    1197
 6 Hertford     1452     7     954  1838     5    1237
 7 Camden        286     0     115   350     2     139
 8 Gates         420     0     254   594     2     371
 9 Warren        968     4     748  1190     2     844
10 Stokes       1612     1     160  2038     5     176
# … with 90 more rows, and 1 more variable:
#   geometry <MULTIPOINT [°]>

plot(st_geometry(nc), border='grey')
plot(st_geometry(nc_pts), pch=16, cex=0.5, add=TRUE)

Casting - POINT

st_cast(nc, "POINT")

Simple feature collection with 2529 features and 7 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
# A tibble: 2,529 × 8
   NAME  BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
   <chr> <dbl> <dbl>   <dbl> <dbl> <dbl>   <dbl>
 1 Ashe   1091     1      10  1364     0      19
 2 Ashe   1091     1      10  1364     0      19
 3 Ashe   1091     1      10  1364     0      19
 4 Ashe   1091     1      10  1364     0      19
 5 Ashe   1091     1      10  1364     0      19
 6 Ashe   1091     1      10  1364     0      19
 7 Ashe   1091     1      10  1364     0      19
 8 Ashe   1091     1      10  1364     0      19
 9 Ashe   1091     1      10  1364     0      19
10 Ashe   1091     1      10  1364     0      19
# … with 2,519 more rows, and 1 more variable:
#   geometry <POINT [°]>

plot(st_geometry(nc), border='grey')
plot(st_geometry(st_cast(nc, "POINT")), pch=16, cex=0.5, add=TRUE)

Casting - LINESTRING

st_cast(nc, "MULTILINESTRING")

Simple feature collection with 100 features and 7 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
# A tibble: 100 × 8
   NAME        BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
   <chr>       <dbl> <dbl>   <dbl> <dbl> <dbl>   <dbl>
 1 Ashe         1091     1      10  1364     0      19
 2 Alleghany     487     0      10   542     3      12
 3 Surry        3188     5     208  3616     6     260
 4 Currituck     508     1     123   830     2     145
 5 Northampton  1421     9    1066  1606     3    1197
 6 Hertford     1452     7     954  1838     5    1237
 7 Camden        286     0     115   350     2     139
 8 Gates         420     0     254   594     2     371
 9 Warren        968     4     748  1190     2     844
10 Stokes       1612     1     160  2038     5     176
# … with 90 more rows, and 1 more variable:
#   geometry <MULTILINESTRING [°]>

st_cast(nc, "MULTILINESTRING") |> st_geometry() |> plot()

Grouping Features

nc_state = st_union(nc)
plot(nc_state)

More Grouping

( nc_cut = nc |>
   mutate(X = st_centroid(nc) |> st_coordinates() |> (\(x) x[,1])()) |>
   mutate(region = cut(X, breaks = 5)) )

Simple feature collection with 100 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
# A tibble: 100 × 10
   NAME        BIR74 SID74 NWBIR74 BIR79 SID79 NWBIR79
 * <chr>       <dbl> <dbl>   <dbl> <dbl> <dbl>   <dbl>
 1 Ashe         1091     1      10  1364     0      19
 2 Alleghany     487     0      10   542     3      12
 3 Surry        3188     5     208  3616     6     260
 4 Currituck     508     1     123   830     2     145
 5 Northampton  1421     9    1066  1606     3    1197
 6 Hertford     1452     7     954  1838     5    1237
 7 Camden        286     0     115   350     2     139
 8 Gates         420     0     254   594     2     371
 9 Warren        968     4     748  1190     2     844
10 Stokes       1612     1     160  2038     5     176
# … with 90 more rows, and 3 more variables:
#   geometry <MULTIPOLYGON [°]>, X <dbl>, region <fct>

ggplot(nc_cut) +
  geom_sf(aes(fill=region))

dplyr and sf

nc_cut |> 
  group_by(region) |> 
  summarize() |> 
  ggplot() + 
    geom_sf(aes(fill=region))

Affine Transfomations

rotate = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)

ctrd = st_centroid(nc_state)
state_rotate = (nc_state) * rotate(-pi/4)
plot(state_rotate, axes=TRUE)

Scaling Size

ctrd = st_centroid(st_geometry(nc))
area = st_area(nc) |> strip_attrs()

nc_rot = nc
st_geometry(nc_rot) = (st_geometry(nc) - ctrd) * 0.75 + ctrd

plot(nc_rot[,"SID79"])

Highway Example

Highways

hwy = read_sf("data/us_interstates.gpkg", quiet = TRUE) |>
  st_transform(st_crs(nc))

ggplot() +
  geom_sf(data=nc) +
  geom_sf(data=hwy, col='red')

NC Interstate Highways

hwy_nc = st_intersection(hwy, nc)

ggplot() + 
  geom_sf(data=nc) +
  geom_sf(data=hwy_nc, col='red')

Counties near the interstate (Projection)

nc_utm  = st_transform(nc, "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs")

ggplot() + 
  geom_sf(data=nc_utm) +
  geom_sf(data=hwy_nc, col='red')

Counties near the interstate (Buffering)

hwy_nc_buffer = hwy_nc |>
  st_transform("+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs") |>
  st_buffer(10000)

ggplot() + 
  geom_sf(data=nc_utm) +
  geom_sf(data=hwy_nc, color='red') +
  geom_sf(data=hwy_nc_buffer, fill='red', alpha=0.3)

Counties near the interstate (Buffering + Union)

hwy_nc_buffer = hwy_nc |>
  st_transform("+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs") |>
  st_buffer(10000) |>
  st_union() |>
  st_sf()
  
ggplot() + 
  geom_sf(data=nc_utm) +
  geom_sf(data=hwy_nc, color='red') +
  geom_sf(data=hwy_nc_buffer, fill='red', alpha=0.3)

Example

How many counties in North Carolina are within 5, 10, 20, or 50 km of an interstate highway?

Gerrymandering Example

NC House Districts - 112th Congress

( nc_house = read_sf("data/nc_districts112.gpkg", quiet = TRUE) |>
    select(ID, DISTRICT) )

Simple feature collection with 13 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32187 ymin: 33.84452 xmax: -75.45998 ymax: 36.58812
Geodetic CRS:  WGS 84
# A tibble: 13 × 3
   ID           DISTRICT                                geom
   <chr>        <chr>                     <MULTIPOLYGON [°]>
 1 037108112001 1        (((-77.32845 35.35031, -77.35398 3…
 2 037108112002 2        (((-78.89928 35.12619, -78.89763 3…
 3 037108112003 3        (((-75.68266 35.23291, -75.68113 3…
 4 037108112004 4        (((-78.77926 35.78568, -78.77947 3…
 5 037108112005 5        (((-79.8968 36.38075, -79.89213 36…
 6 037108112006 6        (((-80.4201 35.68953, -80.41483 35…
 7 037108112007 7        (((-77.59169 34.40907, -77.58699 3…
 8 037108112008 8        (((-78.93373 34.95909, -78.94074 3…
 9 037108112009 9        (((-80.93058 35.18181, -80.9244 35…
10 037108112010 10       (((-81.04032 35.40447, -81.30021 3…
11 037108112011 11       (((-84.00768 35.37262, -84.00831 3…
12 037108112012 12       (((-80.4996 35.6493, -80.51926 35.…
13 037108112013 13       (((-79.90775 36.3818, -79.90694 36…

nc_house = nc_house |>
  st_transform("+proj=utm +zone=17 +datum=NAD83 +units=km +no_defs")
plot(nc_house[,"DISTRICT"], axes=TRUE)

Measuring Compactness - Reock Score

The Reock score is a measure of compactness that is calculated as the the ratio area of a shape to the area of its minimum bounding circle.

circs = nc_house |> 
  lwgeom::st_minimum_bounding_circle()

plot(circs |> filter(DISTRICT == 1) |> st_geometry(), axes=TRUE)
plot(nc_house |> select(DISTRICT) |> filter(DISTRICT == 1), add=TRUE)

ggplot(mapping = aes(fill=DISTRICT)) +
  geom_sf(data=nc_house) +
  geom_sf(data=circs, alpha=0.15) +
  guides(color="none", fill="none")

Calculating Reock

nc_house |>
  mutate(reock = (st_area(nc_house) / st_area(circs)) |> as.numeric()) |>
  ggplot(aes(fill = reock)) +
    geom_sf()

nc_house |>
  mutate(reock = st_area(nc_house) / st_area(circs)) |>
  arrange(reock) |>
  print(n=13)

Simple feature collection with 13 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 196.7724 ymin: 3748.368 xmax: 1002.17 ymax: 4057.317
CRS:           +proj=utm +zone=17 +datum=NAD83 +units=km +no_defs
# A tibble: 13 × 4
   ID           DISTRICT                          geom reock
   <chr>        <chr>              <MULTIPOLYGON [km]>   [1]
 1 037108112012 12       (((545.2987 3945.168, 543.52… 0.116
 2 037108112013 13       (((597.9665 4026.852, 598.03… 0.237
 3 037108112003 3        (((984.0709 3911.853, 984.21… 0.266
 4 037108112002 2        (((691.4141 3889.056, 691.55… 0.303
 5 037108112009 9        (((506.3207 3893.208, 506.88… 0.339
 6 037108112008 8        (((688.6584 3870.456, 688.02… 0.342
 7 037108112011 11       (((226.7522 3918.52, 226.716… 0.344
 8 037108112006 6        (((552.4691 3949.669, 552.94… 0.378
 9 037108112001 1        (((833.677 3918.083, 831.446… 0.378
10 037108112005 5        (((598.9504 4026.746, 599.38… 0.399
11 037108112010 10       (((496.339 3917.9, 472.7444 … 0.411
12 037108112004 4        (((700.7063 3962.453, 700.71… 0.480
13 037108112007 7        (((813.3004 3812.784, 813.74… 0.624

Raster Data (stars)

Example data - Meuse

( meuse_rast = stars::read_stars(
    system.file("external/test.grd", package="raster")
  ) |>
    st_transform(st_crs(meuse_riv))
 )

stars object with 2 dimensions and 1 attribute
attribute(s):
              Min.  1st Qu.   Median    Mean  3rd Qu.
test.grd  138.7071 293.9575 371.9001 425.606 501.0102
              Max. NA's
test.grd  1736.058 6022
dimension(s):
  from  to offset delta              refsys point
x    1  80     NA    NA Amersfoort / RD New    NA
y    1 115     NA    NA Amersfoort / RD New    NA
                      values x/y
x [80x115] 178451,...,181611 [x]
y [80x115] 329530,...,334090 [y]
curvilinear grid

stars class

str(meuse_rast)

List of 1
 $ test.grd: num [1:80, 1:115] NA NA NA NA NA NA NA NA NA NA ...
 - attr(*, "dimensions")=List of 2
  ..$ x:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : num 80
  .. ..$ offset: num NA
  .. ..$ delta : num NA
  .. ..$ refsys:List of 2
  .. .. ..$ input: chr "EPSG:28992"
  .. .. ..$ wkt  : chr "PROJCRS[\"Amersfoort / RD New\",\n    BASEGEOGCRS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            E"| __truncated__
  .. .. ..- attr(*, "class")= chr "crs"
  .. ..$ point : logi NA
  .. ..$ values: num [1:80, 1:115] 178451 178491 178531 178571 178611 ...
  .. ..- attr(*, "class")= chr "dimension"
  ..$ y:List of 7
  .. ..$ from  : num 1
  .. ..$ to    : num 115
  .. ..$ offset: num NA
  .. ..$ delta : num NA
  .. ..$ refsys:List of 2
  .. .. ..$ input: chr "EPSG:28992"
  .. .. ..$ wkt  : chr "PROJCRS[\"Amersfoort / RD New\",\n    BASEGEOGCRS[\"Amersfoort\",\n        DATUM[\"Amersfoort\",\n            E"| __truncated__
  .. .. ..- attr(*, "class")= chr "crs"
  .. ..$ point : logi NA
  .. ..$ values: num [1:80, 1:115] 334090 334090 334090 334090 334090 ...
  .. ..- attr(*, "class")= chr "dimension"
  ..- attr(*, "raster")=List of 4
  .. ..$ affine     : num [1:2] 0 0
  .. ..$ dimensions : chr [1:2] "x" "y"
  .. ..$ curvilinear: logi TRUE
  .. ..$ blocksizes : NULL
  .. ..- attr(*, "class")= chr "stars_raster"
  ..- attr(*, "class")= chr "dimensions"
 - attr(*, "class")= chr "stars"

Plotting

plot(meuse_rast)

plot(
  meuse_riv, 
  col=adjustcolor("lightblue",alpha.f = 0.5), border=NA,
  ylim = c(329500, 333611), axes=TRUE
)
plot(meuse_rast, add=TRUE)

ggplot

ggplot() +
  stars::geom_stars(data=meuse_rast)

ggplot() +
  stars::geom_stars(data=meuse_rast) +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA, alpha=0.5) +
  ylim(329500, 333611)

ggplot() +
  stars::geom_stars(data=meuse_rast) +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA, alpha=0.5) +
  ylim(329500, 333611) +
  scale_fill_gradient(na.value = NA)

Rasters and Projections

meuse_rast_ll = st_transform(meuse_rast, "+proj=longlat +datum=NAD83 +no_defs")

plot(meuse_rast, axes=TRUE)

plot(meuse_rast_ll, axes=TRUE)

meuse_rast

stars object with 2 dimensions and 1 attribute
attribute(s):
              Min.  1st Qu.   Median    Mean  3rd Qu.
test.grd  138.7071 293.9575 371.9001 425.606 501.0102
              Max. NA's
test.grd  1736.058 6022
dimension(s):
  from  to offset delta              refsys point
x    1  80     NA    NA Amersfoort / RD New    NA
y    1 115     NA    NA Amersfoort / RD New    NA
                      values x/y
x [80x115] 178451,...,181611 [x]
y [80x115] 329530,...,334090 [y]
curvilinear grid

meuse_rast_ll

stars object with 2 dimensions and 1 attribute
attribute(s):
              Min.  1st Qu.   Median    Mean  3rd Qu.
test.grd  138.7071 293.9575 371.9001 425.606 501.0102
              Max. NA's
test.grd  1736.058 6022
dimension(s):
  from  to offset delta                       refsys point
x    1  80     NA    NA +proj=longlat +datum=NAD8...    NA
y    1 115     NA    NA +proj=longlat +datum=NAD8...    NA
                        values x/y
x [80x115] 5.72143,...,5.76674 [x]
y [80x115] 50.9557,...,50.9968 [y]
curvilinear grid

Cropping

meuse_rast_riv = meuse_rast[ meuse_riv ]
plot(meuse_rast_riv, axes=TRUE)

Segmentation

meuse_rast |>
  mutate(bins = cut(test.grd, 3) ) |>
  select(bins) |>
  plot()

Polygonization

meuse_rast_poly = meuse_rast |>
  mutate(bins = cut(test.grd, 3) ) |>
  select(bins) |>
  st_as_sf()
plot(meuse_rast_poly)

meuse_rast_poly |>
  group_by(bins) |>
  summarize() |>
  plot()

meuse_rast_poly |>
  group_by(bins) |>
  summarize() |>
  mutate(area = st_area(geometry))

Simple feature collection with 3 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 178431 ymin: 329589.8 xmax: 181590.9 ymax: 333749.8
Projected CRS: Amersfoort / RD New
# A tibble: 3 × 3
  bins                                       geometry   area
* <fct>                                <GEOMETRY [m]>  [m^2]
1 (137,671]          POLYGON ((178551 329829.8, 1785… 4.56e6
2 (671,1.2e+03]      MULTIPOLYGON (((178711 329829.8… 4.74e5
3 (1.2e+03,1.74e+03] MULTIPOLYGON (((179790.9 332189… 5.12e4