skip to Main Content

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

    > library(leaflet)
    > 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:

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

2

Answers


  1. It looks like some of the coordinates in NL have 3 and not 2 coordinates. Row 84 if the first row to show the error. Compare NL$geometry[83] with NL$geometry[84]. The first point for in row 84 is: (-102.3898 22.67024 0).
    I am not sure why there is an extra column in the data.
    Here is a brute force method to separate the differences.

    url = "https://native-land.ca/coordinates/indigenousTerritories.json" 
    NL <- sf::st_read(url)
    
    #dig into the geometry column and determine the number of columns
    columncount <- sapply(1:nrow(NL), function(i) {ncol(NL$geometry[i][[1]][[1]][[1]]) })
    
    #determine the row indexes which has 2 coordinates or 3
    two_col <- which(columncount==2)
    three_col <- which(columncount==3)
    
    #plot the 2 coordinate rows.
    ggplot(NL[two_col,]) +
       geom_sf(aes(geometry = geometry))
    

    I am sure there is a nice sf method to determine the number of coordinate columns but I am not the familiar with the package. As it turns out there 2008 rows with a long&latitude and 51 row with long, lat and other.

    Here is the plot of the 2008 rows:
    enter image description here

    Login or Signup to reply.
  2. 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(ggplot2)
    
    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

    Login or Signup to reply.
Please signup or login to give your own answer.
Back To Top
Search