ateucher / rmapshaper

An R wrapper for the mapshaper javascript library
http://andyteucher.ca/rmapshaper/
Other
200 stars 13 forks source link

ms_simplify: JSON error with mapshaper #111

Closed fm429 closed 3 years ago

fm429 commented 3 years ago

Hi!

I've been trying to speed up some analyses by simplifying an sf data.frame using ms_simplify. The first attempts ended with a few aborted sessions but, thanks to previous issues (#83 and #101), I was able to solve the error by installing nodejs and mapshaper (in cmd.exe as admin).

However, when running the following line of code, I'm getting another error after a few minutes of running time.

UK_simp <- ms_simplify(UK_sel, keep = 0.2, keep_shapes = T, sys = T)

[i] Error: Too-long string in JSON at position 264290734
Run mapshaper -h to view help
Error: Cannot open "C:\...\Temp\RtmpqUySof\file1e0b4782716c3.geojson"; The file doesn't seem to exist.

As I'm not familiar with either mapshaper or JSON, I'm struggling to solve this one. Any suggestion? Thank you in advance!

Filippo

ateucher commented 3 years ago

Hi @fm429 is it possible for you to share the data you are working with so I can try to reproduce this?

fm429 commented 3 years ago

Hi @ateucher, the data I'm working with is from the vector Corine Land Cover 2018 (here's the source link, https://land.copernicus.eu/pan-european/corine-land-cover/clc2018). I clipped it in QGIS to work on a smaller vector file focused on the UK. It's quite a large CSV file and I tried to compress it and upload it here, but unsuccessfully. Is there any other way I can share it? Sorry for the naif question!

ateucher commented 3 years ago

Thanks @fm429 - which format did you download? The geopackage (SQLite) or the ESRI file geodatabase?

fm429 commented 3 years ago

I downloaded the ESRI geodatabase, thanks @ateucher!

ateucher commented 3 years ago

@fm429 I'm unable to reproduce the error - can you please try again with a development version of rmapshaper? You can also avoid using QGIS to select an area by filtering to the UK using a bounding box as you read it into R (I used the geopackage, but should also work with gdb):

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
library(rmapshaper)

gpkg_file <- "~/Desktop/u2018_clc2018_v2020_20u1_geoPackage/DATA/U2018_CLC2018_V2020_20u1.gpkg"
lyr <- "U2018_CLC2018_V2020_20u1"

# Read CRS from file - use SQL to limit to one record and extract the CRS
clc_crs <- st_read(gpkg_file, layer = lyr, query = paste0("SELECT * FROM ", lyr, " LIMIT 1")) %>% 
  st_crs()
#> Reading layer `U2018_CLC2018_V2020_20u1' from data source `/Users/ateucher/Desktop/u2018_clc2018_v2020_20u1_geoPackage/DATA/U2018_CLC2018_V2020_20u1.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 5 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 1916443 ymin: 942165.1 xmax: 1919376 ymax: 943769.9
#> Projected CRS: ETRS89-extended / LAEA Europe

# Create a bbox of the UK in lat/long
uk_bbox_4236 <- st_bbox(c(ymin = 49.674, xmin= -14.015517, ymax= 61.061, xmax = 2.0919117), crs = 4326)

# transform to EPSG:3035 and convert to wkt
uk_wkt <- st_as_text(st_transform(st_as_sfc(uk_bbox_4236), crs = clc_crs))
uk_wkt
#> [1] "POLYGON ((2622823 3232885, 3751375 2981956, 3893574 4241754, 3049221 4437070, 2622823 3232885))"

# Read it in, using the wkt to filter to UK
uk <- st_read(gpkg_file, layer = lyr, wkt_filter = uk_wkt)
#> Reading layer `U2018_CLC2018_V2020_20u1' from data source `/Users/ateucher/Desktop/u2018_clc2018_v2020_20u1_geoPackage/DATA/U2018_CLC2018_V2020_20u1.gpkg' using driver `GPKG'
#> Simple feature collection with 97034 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 2899311 ymin: 2478561 xmax: 4223202 ymax: 4492373
#> Projected CRS: ETRS89-extended / LAEA Europe
summary(uk)
#>    Code_18             Remark             Area_Ha             ID           
#>  Length:97034       Length:97034       Min.   :      0   Length:97034      
#>  Class :character   Class :character   1st Qu.:     37   Class :character  
#>  Mode  :character   Mode  :character   Median :     61   Mode  :character  
#>                                        Mean   :    815                     
#>                                        3rd Qu.:    127                     
#>                                        Max.   :4450136                     
#>            Shape      
#>  MULTIPOLYGON :97034  
#>  epsg:3035    :    0  
#>  +proj=laea...:    0  
#>                       
#>                       
#> 
object.size(uk)
#> 477043848 bytes

# Simplify
uk_simp <- ms_simplify(uk, keep = 0.2, keep_shapes = T, sys = T)

summary(uk_simp)
#>    Code_18             Remark             Area_Ha             ID           
#>  Length:97034       Length:97034       Min.   :      0   Length:97034      
#>  Class :character   Class :character   1st Qu.:     37   Class :character  
#>  Mode  :character   Mode  :character   Median :     61   Mode  :character  
#>                                        Mean   :    815                     
#>                                        3rd Qu.:    127                     
#>                                        Max.   :4450136                     
#>            Shape      
#>  MULTIPOLYGON :97034  
#>  epsg:3035    :    0  
#>  +proj=laea...:    0  
#>                       
#>                       
#> 
object.size(uk_simp)
#> 169216552 bytes

plot(uk_simp[, "Code_18"])

Created on 2021-05-20 by the reprex package (v2.0.0)

ateucher commented 3 years ago

I forgot to mention: you can install the dev version of mapshaper with:

# install.packages("remotes")
remotes::install_github("ateucher/rmapshaper", ref = "fix-erase-bug")
fm429 commented 3 years ago

This is brilliant @ateucher! Thank you so much. I've used the resulting file, uk_simp, in a function where I'm retrieving the type of landscape (Code_18) based on a series of coordinates and seems to work fine using st_join(x, uk_simp, join = st_within)). I got an NA value but I'll have a better look when I have time. Thanks again.

As for the dev version of mapshaper, I tried to install it but I got this error: Error in utils::download.file(url, path, method = method, quiet = quiet, : cannot open URL 'https://api.github.com/repos/ateucher/rmapshaper/tarball/fix-erase-bug'.

ateucher commented 3 years ago

Great to hear @fm429. Sorry, I deleted the fix-erase-bug branch. You can now install the dev version with just remotes::install_github("ateucher/rmapshaper")