Read NativeLand JSON and plot in ggplot

I would like to add my own features to the map of Native Land:

    > library(leaflet)
    > library(sf)

    > url = "https://native-land.ca/coordinates/indigenousTerritories.json" 
    > NL <- sf::st_read(url)
    > ggplot(NL) + geom_sf()

But I get the follow error:

Error in `geom_sf()`:
! Problem while converting geom to grob.
ℹ Error occurred in the 1st layer.

When I use this site to validate the GeoJSON, it indicates:

Line 1: Polygons and MultiPolygons should follow the right-hand rule

Any insight to use this dataset is welcome !

Previous efforts:

Not sure why, but this works:

    > readLines(url) %>%
        leaflet() %>% 
        addProviderTiles("Esri.WorldImagery") %>% 
        addGeoJSON(native_land_json)

Solution

Just to add to the previous answer. The excluded section of that ggplot error message, along with the backtrace, would give a few valuable hints:

Caused by error:
! number of columns of matrices must match (see arg 84)
...
 24.                           └─sf:::st_as_grob.sfc_MULTIPOLYGON(...)
 25.                             ├─base::do.call(rbind, x)

So what could cause rbind to fail with that error? When checking the NL object, it reports Dimension: XY, XYZ and z_range: zmin: 0 zmax: 0, first being an indicator that not all geometries share the same dimension:

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
url = "https://native-land.ca/coordinates/indigenousTerritories.json" 
NL <- read_sf(url)

NL |> 
  print(n = 3)
#> Simple feature collection with 2059 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY, XYZ
#> Bounding box:  xmin: -189.8148 ymin: -56.02359 xmax: 182.1078 ymax: 83.77334
#> z_range:       zmin: 0 zmax: 0
#> Geodetic CRS:  WGS 84
#> # A tibble: 2,059 × 7
#>   id               Name     ID Slug  description color                  geometry
#>   <chr>            <chr> <int> <chr> <chr>       <chr>        <MULTIPOLYGON [°]>
#> 1 5ca71dcffa94678… Tari… 39506 tari… https://na… #0e1… (((152.6523 -25.15547, 1…
#> 2 4777473e426a4b6… Djiru 39497 djiru https://na… #0e1… (((146.1055 -17.81829, 1…
#> 3 a3f2e55a99b0523… Girr… 39496 girr… https://na… #0e1… (((146.0295 -18.08393, 1…
#> # ℹ 2,056 more rows

And we can see that the number of columns in coordinate matrices are indeed not matching:

st_geometry(NL[83,]) |> st_coordinates() |> head(n=2)
#>             X         Y L1 L2 L3
#> [1,] 124.8816 -18.47961  1  1  1
#> [2,] 124.7498 -18.67747  1  1  1
st_geometry(NL[84,]) |> st_coordinates() |> head(n=2)
#>              X        Y Z L1 L2 L3
#> [1,] -102.3898 22.67024 0  1  1  1
#> [2,] -102.1236 22.36702 0  1  1  1

To drop Z from all geometries in that set, we can use st_zm(), note the Dimension: XY; resulting sf object is also acceptable for ggplot:

NL |> 
  st_zm(drop = TRUE, what = "ZM") |>
  print(n = 3)
#> Simple feature collection with 2059 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -189.8148 ymin: -56.02359 xmax: 182.1078 ymax: 83.77334
#> Geodetic CRS:  WGS 84
#> # A tibble: 2,059 × 7
#>   id               Name     ID Slug  description color                  geometry
#> * <chr>            <chr> <int> <chr> <chr>       <chr>        <MULTIPOLYGON [°]>
#> 1 5ca71dcffa94678… Tari… 39506 tari… https://na… #0e1… (((152.6523 -25.15547, 1…
#> 2 4777473e426a4b6… Djiru 39497 djiru https://na… #0e1… (((146.1055 -17.81829, 1…
#> 3 a3f2e55a99b0523… Girr… 39496 girr… https://na… #0e1… (((146.0295 -18.08393, 1…
#> # ℹ 2,056 more rows

One option to detect/subset individual geometries with Z dimension is through checking class inheritance (XY vs XYZ) :

is_xyz <- function(sfc) {
  sapply(sfc, inherits, "XYZ")
}

# XYZ subset
NL[ is_xyz(st_geometry(NL)), ] |> print(n = 3)
#> Simple feature collection with 51 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XYZ
#> Bounding box:  xmin: -168.5537 ymin: -33.0908 xmax: -39.98129 ymax: 68.51203
#> z_range:       zmin: 0 zmax: 0
#> Geodetic CRS:  WGS 84
#> # A tibble: 51 × 7
#>   id               Name     ID Slug  description color                  geometry
#>   <chr>            <chr> <int> <chr> <chr>       <chr>        <MULTIPOLYGON [°]>
#> 1 01287d68c58c508… Huac… 36648 huac… https://na… #B9E… Z (((-102.3898 22.67024 …
#> 2 188c197984cc05f… Tecu… 36625 tecu… https://na… #D26… Z (((-103.3731 20.93835 …
#> 3 b1801173f6356d9… jñat… 36590 maza… https://na… #68F… Z (((-100.3474 19.74578 …
#> # ℹ 48 more rows

# XY  subset
NL[!is_xyz(st_geometry(NL)), ] |> print(n = 3)
#> Simple feature collection with 2008 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -189.8148 ymin: -56.02359 xmax: 182.1078 ymax: 83.77334
#> Geodetic CRS:  WGS 84
#> # A tibble: 2,008 × 7
#>   id               Name     ID Slug  description color                  geometry
#>   <chr>            <chr> <int> <chr> <chr>       <chr>        <MULTIPOLYGON [°]>
#> 1 5ca71dcffa94678… Tari… 39506 tari… https://na… #0e1… (((152.6523 -25.15547, 1…
#> 2 4777473e426a4b6… Djiru 39497 djiru https://na… #0e1… (((146.1055 -17.81829, 1…
#> 3 a3f2e55a99b0523… Girr… 39496 girr… https://na… #0e1… (((146.0295 -18.08393, 1…
#> # ℹ 2,005 more rows

^{Created on 2024-01-22 with reprex v2.0.2}