R: Recognizing Geographical Coordinates within JSON

I am working with the R programming language.

My Problem: I am trying to recreate this map (e.g. in leaflet) here myself : https://crtc.gc.ca/cartovista/LTEOverTheYearsYE2019_EN/index.html (e.g. for 2019).

Here is what I have tried so far:

Step 1 : First, I downloaded the JSON for 2019:

library(jsonlite)
library(httr)


response <- GET("https://crtc.gc.ca/cartovista/LTEOverTheYearsYE2019_EN/map/LTE_YE2019.json")

content <- content(response, as = "text")

data <- fromJSON(content)

Step 2: Next, I tried to extract the longitude/latitude for all the polygons in this file and convert them into the correct Coordinate Reference System:

library(sp)
library(rgdal)

polygons <- lapply(LTE_YE2019$f$g$c[[1]], function(x) {
  Polygon(matrix(x, ncol = 2, byrow = TRUE))
})

sp_polygons <- SpatialPolygons(list(Polygons(polygons, ID = 1)))


spdf <- SpatialPolygonsDataFrame(sp_polygons, data.frame(id = 1, row.names = 1))

proj4string(spdf) <- CRS(LTE_YE2019$proj)

The file now looks something like this:

> str(spdf)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame':	1 obs. of  1 variable:
  .. ..$ id: num 1
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 606

....

  .. .. .. .. .. [list output truncated]
  .. .. .. ..@ plotOrder: int [1:606] 321 322 324 329 330 325 372 290 267 373 ...
  .. .. .. ..@ labpt    : num [1:2] 994419 305033
  .. .. .. ..@ ID       : chr "1"
  .. .. .. ..@ area     : num 1.17e+11
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] -2319216 -307601 2981680 3014094
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=lcc +lat_0=49 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"
  .. .. ..$ comment: chr "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"Unknown based on GRS80 ellipsoid\",\n      "| __truncated__
  ..$ comment: chr "FALSE"

My Question: When I try to plot this, I get the following error:

library(leaflet)

m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = spdf)

#Error in rgeos::createPolygonsComment(pgons) : 
  #rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 2

Step 3: To solve this problem, I am trying to following approaches :

# Approach 1:

require(devtools)
install_github("eblondel/cleangeo")
require(cleangeo)

    report <- clgeo_CollectionReport(spdf)
    summary <- clgeo_SummaryReport(report)
    issues <- report[report$valid == FALSE,]
    
    spdf.clean <- clgeo_Clean(spdf)


# Approach 2 

library(rmapshaper)

spdf_clean <- ms_simplify(spdf)

m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = spdf_clean)

# Approach 3

library(sf)
library(leaflet)


sf_polygons <- st_as_sf(spdf)
valid <- st_is_valid(sf_polygons)
sf_polygons_valid <- st_make_valid(sf_polygons)


m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = sf_polygons_valid)

All 3 approaches seem to be taking a long time to run - in the meantime, I am not sure if I am approaching this problem correctly.

Can someone please tell me if I am doing this correctly?

How can I resolve this error to recreate this map in R/leaflet?

Thanks!

Note: Perhaps I have significantly overcomplicated this question. In the end, my main goal is to replicate these geographical maps as interactive leaflet maps in R - perhaps there might be an easier and more direct way to do this? I am open any suggestions!

Full credit for JSON parsing goes to Spacedman, his answer - https://gis.stackexchange.com/a/461914/213352 - to the related question basically solved the current one. The following is just a purrr-heavy attempt to pack it all together with some questionable map hole handling thrown in -- some of those polygons are holes and some are islands in holes, probably with matching ring directions. Here holes were handled through rather crude combination of predicates, unions and removals.

library(jsonlite)
library(purrr)
library(dplyr)
library(sf)
library(sfheaders)
library(leaflet)

vect_to_sfg <- function(coord_vect){
  # input vector of differential coordinates, 
  # [1] & [2] are refernces for x & y:
  # int [1:18242] 2016758 388079 -1045 -494 -1046 -494 ...
  
  # convert to 2 column matrix
  #          [,1]   [,2]
  # [1,] 2016758 388079
  # [2,]   -1045   -494
  # [3,]   -1046   -494  
  
  # apply cumsum to columns to transform differential values to 
  # absolute (x,y) coordinates
  #         [,1]   [,2]
  # [1,] 2016758 388079
  # [2,] 2015713 387585
  # [3,] 2014667 387091
  
  # convert to sfg, a single polygon geomtery
  matrix(coord_vect, ncol = 2, byrow = TRUE) %>%  
    apply(2, cumsum) %>% 
    sfheaders::sfg_polygon()
}

# list of vectors ( -> list of polygons ) -> sf geometry column
vect_list_to_sfc <- function (coord_vect_list) {
  map(coord_vect_list, vect_to_sfg) %>% 
    st_sfc()
}

# crude approach to cut out holes
make_holes <- function(sfc){
  p_union      <- st_union(sfc)
  hole_mask    <- st_contains_properly(p_union, sfc, sparse = FALSE)
  hole_union   <- st_union(sfc[hole_mask])
  island_mask  <- st_contains_properly(hole_union, sfc, sparse = FALSE)
  island_union <- st_union(sfc[island_mask])
  st_difference(p_union, hole_union) %>% st_union(island_union)
}

#### fetch and process

base_url <- "https://crtc.gc.ca/cartovista/LTEOverTheYearsYE2019_EN/map/"
jsons <- paste0("LTE_YE", 2014:2019, ".json")

# request and parse jsons, returns list of (proj string; geometry)
j_resp <- set_names(jsons) %>% 
  map(\(json_f) fromJSON(paste0(base_url, json_f))) %>% 
  map(\(j) list(proj = pluck(j, "proj"),
                geom = pluck(j, "f", "g", "c", 1)))

# sample from j_resp:
glimpse(j_resp)
#> List of 6
#>  $ LTE_YE2014.json:List of 2
#>   ..$ proj: chr "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs "
#>   ..$ geom:List of 265
#>   .. ..$ : int [1:18242] 2016758 388079 -1045 -494 -1046 -494 -1045 -493 -1045 -493 ...
#>   .. ..$ : int [1:2182] 1472249 -197163 -1189 -360 -1188 -360 -1189 -360 -1188 -360 ...
#>   .. .. [list output truncated]

# extract proj strimgs
projs <- map(j_resp, "proj")

# extract geometry vector lists, 
# convert to sf geometry columns, 
# handle holes (overlapping geometries),
# set proj for every sf object,
# transform to WGS84 for leaflet

# result is a list of sf objects, one multipolygon per each input json

lte_sf_lst <- j_resp %>% 
  map("geom") %>% 
  map(vect_list_to_sfc) %>% 
  map(make_holes) %>% 
  map2(projs, \(geom, proj) st_sf(geometry = geom, crs = proj)) %>% 
  map(st_transform, crs = "WGS84")

# helper for reduce2, add new polygon layer to existing leaflet object
addPolygons_reduce <-  function(map, polygon, name){
  addPolygons(map, data = polygon, group = name, weight = 0, 
              fillColor = "#3b528b", fillOpacity = .2)
}

# inital Leaflet object, no layers, serves as .init argument for reduce2
lte_leaflet <- leaflet() %>% 
  addProviderTiles(providers$CartoDB.Positron)

# use reduce2 to add all collected layers to Leaflet,
# iterate through lte_sf_lst list and call addPolygons_reduce() on 
# each object, initial map being lte_leaflet
lte_leaflet <- reduce2(lte_sf_lst, 
                       names(lte_sf_lst), 
                       addPolygons_reduce, 
                       .init = lte_leaflet) %>% 
  # add layer contorls
  addLayersControl(overlayGroups = names(lte_sf_lst), 
                   options = layersControlOptions(collapsed = FALSE))
lte_leaflet

<!-- -->

^{Created on 2023-06-30 with reprex v2.0.2}