Add utility functions for splitting and weighting LINESTRING and POLYGON input geometry

[elipousson]

Following up on our call last week, here is a first rough draft at a utility function with a reprex to split an input sf object and assign a weight based on length or physical area. This could be modified to make a split resource by Census tract that could use block-level population weighting to determine the weight assigned to the resource in one tract vs. another.

The reprex uses my mapbaltimore package just for convenience so I can rework with more generic input data if that is helpful.

library(tidyverse)
library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(rlang)
#> 
#> Attaching package: 'rlang'
#> The following objects are masked from 'package:purrr':
#> 
#>     %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
#>     flatten_raw, invoke, splice

split_resource_by_intersection <- function(
    resource,
    by = NULL,
    weight = c("length", "area"),
    ...) {
  stopifnot(
    inherits(resource, "sf"),
    inherits(by, c("sfc", "sf")),
    any(sf::st_is(by, c("POLYGON", "MULTIPOLYGON")))
  )
  
  crs <- sf::st_crs(resource)
  
  if (sf::st_crs(resource) != sf::st_crs(by)) {
    by <- sf::st_transform(by, crs = crs)
  }
  
  weight <- arg_match0(weight, c("length", "area"))
  
  if (weight == "length") {
    stopifnot(
      any(sf::st_is(resource, c("LINESTRING", "MULTILINESTRING")))
    )
    weight_fn <- sf::st_length
  } else if (weight == "area") {
    stopifnot(
      any(sf::st_is(resource, c("POLYGON", "MULTIPOLYGON")))
    )
    weight_fn <- sf::st_area
  }
  
  resource_intersection <- sf::st_intersection(
    x = resource,
    y = by,
    ...
  )
  
  stopifnot(
    !has_name(resource, weight)
  )
  
  resource_intersection <- resource_intersection |> 
    dplyr::mutate(
      "{weight}" := suppressWarnings(weight_fn(geometry))
    )
  
  
  resource_intersection
}

length_split <- split_resource_by_intersection(
  mapbaltimore::mta_bus_lines,
  mapbaltimore::baltimore_tracts,
  "length"
)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

length_split
#> Simple feature collection with 1150 features and 15 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: 424876.3 ymin: 170183.3 xmax: 440585.1 ymax: 189421
#> Projected CRS: NAD83(HARN) / Maryland
#> First 10 features:
#>                           route_name                route_type    route_number
#> 46       WOODBERRY - CANTON CROSSING MTA Local Bus - LocalLink              21
#> 47               MONDAWMIN - BAYVIEW MTA Local Bus - LocalLink              22
#> 57                 TOWSON - DOWNTOWN MTA Local Bus - LocalLink              51
#> 86     FORT MCHENRY - SINAI HOSPITAL MTA Local Bus - LocalLink              94
#> 87            DOWNTOWN - ROLAND PARK MTA Local Bus - LocalLink              95
#> 98       CURTIS BAY - HOPKINS/MORGAN  MTA Local Bus - CityLink CityLink SILVER
#> 81         PENN-NORTH - MILFORD MILL MTA Local Bus - LocalLink              85
#> 47.1             MONDAWMIN - BAYVIEW MTA Local Bus - LocalLink              22
#> 58   GREENMOUNT NORTH - STELLA MARIS MTA Local Bus - LocalLink              52
#> 97     DOWNTOWN - TOWSON/LUTHERVILLE  MTA Local Bus - CityLink    CityLink RED
#>      route_abb frequent school tractce       geoid    name             namelsad
#> 46          21    FALSE  FALSE  120202 24510120202 1202.02 Census Tract 1202.02
#> 47          22     TRUE  FALSE  120202 24510120202 1202.02 Census Tract 1202.02
#> 57          51    FALSE  FALSE  120202 24510120202 1202.02 Census Tract 1202.02
#> 86          94    FALSE  FALSE  120202 24510120202 1202.02 Census Tract 1202.02
#> 87          95    FALSE  FALSE  120202 24510120202 1202.02 Census Tract 1202.02
#> 98          SV     TRUE  FALSE  120202 24510120202 1202.02 Census Tract 1202.02
#> 81          85     TRUE  FALSE  272005 24510272005 2720.05 Census Tract 2720.05
#> 47.1        22     TRUE  FALSE  120201 24510120201 1202.01 Census Tract 1202.01
#> 58          52    FALSE  FALSE  120201 24510120201 1202.01 Census Tract 1202.01
#> 97          RD     TRUE  FALSE  120201 24510120201 1202.01 Census Tract 1202.01
#>       aland awater    intptlat     intptlon                       geometry
#> 46   981794      0 +39.3286922 -076.6194426 LINESTRING (432337.1 184235...
#> 47   981794      0 +39.3286922 -076.6194426 MULTILINESTRING ((432414.2 ...
#> 57   981794      0 +39.3286922 -076.6194426 MULTILINESTRING ((432909.8 ...
#> 86   981794      0 +39.3286922 -076.6194426 LINESTRING (432310.5 184276...
#> 87   981794      0 +39.3286922 -076.6194426 MULTILINESTRING ((433049.8 ...
#> 98   981794      0 +39.3286922 -076.6194426 MULTILINESTRING ((432997.2 ...
#> 81   927829      0 +39.3686872 -076.7010065 MULTILINESTRING ((425872.8 ...
#> 47.1 348020      0 +39.3279849 -076.6118803 MULTILINESTRING ((433676.1 ...
#> 58   348020      0 +39.3279849 -076.6118803 MULTILINESTRING ((433673.5 ...
#> 97   348020      0 +39.3279849 -076.6118803 MULTILINESTRING ((433673.5 ...
#>              length
#> 46     49.17965 [m]
#> 47    689.82488 [m]
#> 57   2017.86601 [m]
#> 86     49.17965 [m]
#> 87   2471.48348 [m]
#> 98   2563.20195 [m]
#> 81   1370.39830 [m]
#> 47.1  874.10713 [m]
#> 58    135.89796 [m]
#> 97    135.89796 [m]

area_split <- split_resource_by_intersection(
  mapbaltimore::neighborhoods,
  mapbaltimore::baltimore_tracts,
  "area"
)
#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

area_split
#> Simple feature collection with 1161 features and 14 fields
#> Geometry type: GEOMETRY
#> Dimension:     XY
#> Bounding box:  xmin: 424869.7 ymin: 170004.4 xmax: 440590.2 ymax: 189405.1
#> Projected CRS: NAD83(HARN) / Maryland
#> # A tibble: 1,161 × 15
#>    name  type  acres osm_id wikidata tractce geoid name.1 namelsad  aland awater
#>  * <chr> <chr> <dbl> <chr>  <chr>    <chr>   <chr> <chr>  <chr>     <dbl>  <dbl>
#>  1 Abell Resi…  46.7 12740… Q4666688 120202  2451… 1202.… Census … 981794      0
#>  2 Char… Resi… 169.  12740… Q5083136 120202  2451… 1202.… Census … 981794      0
#>  3 Guil… Resi… 354.  12740… Q5615899 120202  2451… 1202.… Census … 981794      0
#>  4 John… Resi… 164.  <NA>   <NA>     120202  2451… 1202.… Census … 981794      0
#>  5 Oake… Resi…  40.0 12740… Q7073836 120202  2451… 1202.… Census … 981794      0
#>  6 Remi… Resi… 157.  12755… Q7311815 120202  2451… 1202.… Census … 981794      0
#>  7 Wyma… Resi…  70.6 12744… Q8039970 120202  2451… 1202.… Census … 981794      0
#>  8 Ches… Resi… 547.  12800… Q5094038 272005  2451… 2720.… Census … 927829      0
#>  9 Cros… Resi… 379.  12800… <NA>     272005  2451… 2720.… Census … 927829      0
#> 10 Fall… Resi… 189.  12802… Q5432552 272005  2451… 2720.… Census … 927829      0
#> # ℹ 1,151 more rows
#> # ℹ 4 more variables: intptlat <chr>, intptlon <chr>, geometry <GEOMETRY [m]>,
#> #   area [m^2]

^{Created on 2024-09-24 with reprex v2.1.1}

[elipousson]

I worked up this draft into a more developed version in this feature branch: https://github.com/elipousson/sedtR/blob/weight-resource-feat/R/weight_resource.R

Feedback welcome!

[elipousson]

Hey @Deckart2! The linked branch is where the weight functions are located. I made some further updates to this branch based on the discussion for the multi-format PR including setting a preference for sf::st_point_on_surface and adding some tests.

I can open this as a PR after the multi-format support PR and the visualization PR is finished (or merge it into the multi-format branch if you want to do it all in one).

[Deckart2]

Ah thanks @elipousson! I think keeping them as separate makes sense for now.