tmap_options(check.and.fix = TRUE) doesn't fix invalid geometries unless sf::sf_use_s2(FALSE) is used

[neon-ninja]

Hi,

After updating to the latest version of tmap, some code that was written for an earlier version of tmap no longer works. The issue appears to be caused by invalid geometries. The error is thrown even with tmap_options(check.and.fix = TRUE). For the purposes of a reprex, I've filtered the invalid geometries out of the larger dataset. Here's a reprex:

With invalid_geometry.rds.zip in your working directory:

library(sf)
#> Warning: package 'sf' was built under R version 4.0.5
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

invalid_geometry = readRDS("invalid_geometry.rds.zip")
st_is_valid(invalid_geometry, reason = TRUE)
#> [1] "Edge 9 crosses edge 26"  "Edge 0 crosses edge 2"  
#> [3] "Edge 18 crosses edge 20" "Edge 0 crosses edge 15" 
#> [5] "Edge 3 crosses edge 5"

library(remotes)
#> Warning: package 'remotes' was built under R version 4.0.5
install_github("r-tmap/tmaptools")
#> Skipping install of 'tmaptools' from a github remote, the SHA1 (0c8b0b1c) has not changed since last install.
#>   Use `force = TRUE` to force installation
install_github("r-tmap/tmap")
#> Skipping install of 'tmap' from a github remote, the SHA1 (2e737efa) has not changed since last install.
#>   Use `force = TRUE` to force installation

library(tmap)

tmap_options(check.and.fix = TRUE)

tm_shape(invalid_geometry) +
  tm_fill()
#> Warning: The shape invalid_geometry is invalid. See sf::st_is_valid
#> Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
#>   Evaluation error: Found 5 features with invalid spherical geometry.
#> [1] Loop 0 is not valid: Edge 9 crosses edge 26
#> [2] Loop 0 is not valid: Edge 0 crosses edge 2
#> [3] Loop 0 is not valid: Edge 18 crosses edge 20
#> [4] Loop 0 is not valid: Edge 0 crosses edge 15
#> [5] Loop 0 is not valid: Edge 3 crosses edge 5.
#> Error in co[, 1]: incorrect number of dimensions

^{Created on 2021-10-12 by the reprex package (v2.0.1)}

I also note that st_make_valid doesn't seem to help:

st_is_valid(st_make_valid(invalid_geometry))
[1] FALSE FALSE FALSE FALSE FALSE

cleangeo::clgeo_Clean(as_Spatial(invalid_geometry))

seems to fix it

[Nowosad]

Have you tried to set st_use_s2(FALSE)?

[neon-ninja]

Doesn't look like my version of sf (1.0-2) has that

> sf::st_use_s2(FALSE)
Error: 'st_use_s2' is not an exported object from 'namespace:sf'

[Nowosad]

Sorry - I meant sf_use_s2(FALSE)

[neon-ninja]

That seems to help, it changes it from an error to a warning:

> sf::sf_use_s2(FALSE)
Spherical geometry (s2) switched off
> tm_shape(invalid_geometry) +
+   tm_fill()
Warning message:
The shape invalid_geometry is invalid. See sf::st_is_valid

[mtennekes]

This is mostly beyond tmap, because tmap expects valid geometries. The check.and.fix option uses st_make_valid under the hood.

The fact that st_make_valid(invalid_geometry) doesn't work when sf_use_s2(TRUE) is something that should be reported in the sf repo.

[neon-ninja]

Looks like someone already raised this issue over in the sf repo - r-spatial/sf#1732. Can I suggest tmap falling back to cleangeo::clgeo_Clean(as_Spatial(invalid_geometry)) if st_make_valid fails? Alternatively setting sf::sf_use_s2(FALSE) until st_make_valid is fixed.

[mtennekes]

Thanks for the suggestions.

The first suggestion will be a no, because I don't want to include extra dependencies (cleangeo and sp)
We will think about the second suggestion.

@Nowosad what is your opinion?

[Nowosad]

My current suggestion would be to update the warning or move to the error message (e.g., Warning: The shape invalid_geometry is invalid. Try to apply sf::st_is_valid() or use sf::sf_use_s2(FALSE) before making the map.).

There is also another thing to consider related to s2 - how to plot lines in tmap when sf::sf_use_s2(TRUE)? Should they be straight lines or the shortest paths (see r-spatial/sf#1780 (comment))?

[eblondel]

My 2 cents here, discovering this issue.
if you deal with simple features (sf) that as they are called, are managed or come from tools based on OGC specifications (and maybe how those were interpreted) including a specific frame for assessing the geometry validity; which is the case for most of GIS software (proprietary or from the OsGeo) then the best is to use sf::sf_use_s2(FALSE); On the contrary, which is now the default behavior of sf (we may wonder why), you will rely on the geometry validity rules that come from s2 library.
If you grab data from OGC services in which data is considered as valid, grabing them with sf::sf_use_s2(TRUE) you may end up with non-valid geometries; By sf::sf_use_s2(FALSE), you will rely on the same assertions that the other GIS tools.

[pennybeaver]

Hi Martijn or Jakub

I am currently transferring a thesis to publication, my thesis used tmap for all spatial data analysed which as submitted late 2021. So here I am recreating some of the maps and I am getting the same error, I realise tmap will not be supported or updated. I have tried the above fixes with no change in the error message. I think it is in relation to the readtopo dimension but any suggestions you have will be greatly appreciated.

CODE

library(sf)
library(tmap)
library(tmaptools)
data(World)
tmap_mode("plot")

#Breeding populations locations
coordinates(sites) <- ~ lon + lat
siteMT <- sites[sites$site =="Muttonbird Island",]
siteMG <- sites[sites$site =="Montague Island",]

#obtain bathymetry data for map
bathywedge <- readtopo("etopo2", xylim = extent(130, 180, -60, -25))

bathywedge2<-rasterToContour(bathywedge)
class(bathywedge2)

tmap::tm_shape(bathywedge2)+tm_iso()+
tmap::tm_shape(World)+
tm_polygons(col = "black")+
tm_text("name", bg.color = "black", bg.alpha = 0.5, remove.overlap = T, size.lowerbound = T, scale = 1, size = 0.7, colorNULL = TRUE) +
tmap::tm_shape(islandname)+tm_text("name", size = 0.7, col = "white")+
tm_grid(col ="gray80", alpha = 0.3)+
tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)+
tm_shape(siteMT)+tm_dots(size = 0.5, shape = 16, col = "red")+
tm_shape(siteMG)+tm_dots(size = 0.5, shape = 16, col = "red")+
tm_shape(siteGI)+tm_dots(size = 0.5, shape = 16, col = "orange")+
tm_shape(Heron)+tm_dots(size = 0.5, shape = 16, col = "orange")+
tm_shape(LHI)+tm_dots(size = 0.5, shape = 16, col = "orange")+
tm_shape(NC)+tm_dots(size = 0.5, shape = 16, col = "orange")

Error in cls[2, ] : incorrect number of dimensions
In addition: Warning messages:
1: Currect projection of shape siteMT unknown. Long-lat (WGS84) is assumed.
2: Currect projection of shape siteMG unknown. Long-lat (WGS84) is assumed.
3: Currect projection of shape siteGI unknown. Long-lat (WGS84) is assumed.
4: Currect projection of shape Heron unknown. Long-lat (WGS84) is assumed.
5: Currect projection of shape LHI unknown. Long-lat (WGS84) is assumed.
6: Currect projection of shape NC unknown. Long-lat (WGS84) is assumed.

[mtennekes]

Did you try sf::sf_use_s2(FALSE) @pennybeaver ?
If it still doesn't work a minimal working example that I can reproduce would help.

Planning to submit the last 3.x version this week, so if possible/needed I can fix this.

[pennybeaver]

Did you try sf::sf_use_s2(FALSE) @pennybeaver ? If it still doesn't work a minimal working example that I can reproduce would help.

Planning to submit the last 3.x version this week, so if possible/needed I can fix this.

Yes I have tried sf::sf_use_s2(FALSE) see below minimal working example. When I break it down it is the readtopo code that is the issue.

library(tmap)
library(tmaptools)
data(World)
tmap_mode("plot")
library(sf)
sf::sf_use_s2(FALSE)

bathywedge <- readtopo("etopo2", xylim = extent(130, 180, -60, -25))
class(bathywedge)

bathywedge2<-rasterToContour(bathywedge)
class(bathywedge2)

tmap::tm_shape(bathywedge2)+tm_iso()+
tmap::tm_shape(World)+
tm_polygons(col = "black")+
tm_text("name", bg.color = "black", bg.alpha = 0.5, remove.overlap = T, size.lowerbound = T, scale = 1, size = 0.7, colorNULL = TRUE) +
tm_grid(col ="gray80", alpha = 0.3)+
tm_xlab("Longitude", size = 0.8)+tm_ylab("Latitude", size = 0.8)

[Nowosad]

@pennybeaver I still cannot reproduce the code. Could you maybe save the bathywedge2 object to a file (e.g., using saveRDS, compress it to a zip file, and then attach it here?

[pennybeaver]

[https://res-h3.public.cdn.office.net/assets/mail/file-icon/png/zip_16x16.png]Tmap data files.RData.zip<https://universitytasmania-my.sharepoint.com/:u:/g/personal/penny_beaver_utas_edu_au/ERzHtIwGpC1Frs2vUkAX6fYBdSWJXigRblhXzwv1LJS5Mw> HI Jakub Attached is the zip file. Thanks Penny

…

________________________________ From: Jakub Nowosad ***@***.***> Sent: Tuesday, 5 September 2023 8:09 PM To: r-tmap/tmap ***@***.***> Cc: Penny Beaver ***@***.***>; Mention ***@***.***> Subject: Re: [r-tmap/tmap] Error in co[, 1] : incorrect number of dimensions (#606) @pennybeaver<https://github.com/pennybeaver> I still cannot reproduce the code. Could you maybe save the bathywedge2 object to a file (e.g., using saveRDS, compress it to a zip file, and then attach it here? — Reply to this email directly, view it on GitHub<#606 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AFIRKD7E5F7PY7GXDAKKJIDXY326TANCNFSM5FZG4R5A>. You are receiving this because you were mentioned.Message ID: ***@***.***> This email is confidential, and is for the intended recipient only. Access, disclosure, copying, distribution, or reliance on any of it by anyone outside the intended recipient organisation is prohibited and may be a criminal offence. Please delete if obtained in error and email confirmation to the sender. The views expressed in this email are not necessarily the views of the University of Tasmania, unless clearly intended otherwise.

[pennybeaver]

Tmap data files.RData.zip

zip file attached :)

[Nowosad]

@pennybeaver are you using the CRAN version of tmap or the GitHub one? The second one seems to be working (remotes::install_github("r-tmap/tmap")):

# remotes::install_github("r-tmap/tmap")
library(tmap)
library(tmaptools)
data(World)
tmap_mode("plot")
#> tmap mode set to plotting
library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

load("Tmap data files.RData")

tmap::tm_shape(bathywedge2) + tm_iso() +
  tmap::tm_shape(World) +
  tm_polygons(col = "black") +
  tm_text(
    "name",
    bg.color = "black",
    bg.alpha = 0.5,
    remove.overlap = T,
    size.lowerbound = T,
    scale = 1,
    size = 0.7,
    colorNULL = TRUE
  ) +
  tm_grid(col = "gray80", alpha = 0.3) +
  tm_xlab("Longitude", size = 0.8) + tm_ylab("Latitude", size = 0.8)

[pennybeaver]

I have reinstalled twice once on an external cloud RStudio site (organisational site) and at the same time on the pc hard drive. Both had the error message originally but not the hard drive after I reinstalled tmap using the below code.

remotes::install_github("r-tmap/tmap")

Thank you very much for your help, it is appreciated.

[neon-ninja]

This is still an issue, the error is just slightly different now:

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

invalid_geometry = readRDS("invalid_geometry.rds.zip")
st_is_valid(invalid_geometry, reason = TRUE)
#> [1] "Edge 9 crosses edge 26"  "Edge 0 crosses edge 2"  
#> [3] "Edge 18 crosses edge 20" "Edge 0 crosses edge 15" 
#> [5] "Edge 3 crosses edge 5"

library(remotes)
install_github("r-tmap/tmaptools")
#> Skipping install of 'tmaptools' from a github remote, the SHA1 (0c8b0b1c) has not changed since last install.
#>   Use `force = TRUE` to force installation
install_github("r-tmap/tmap")
#> Skipping install of 'tmap' from a github remote, the SHA1 (40106020) has not changed since last install.
#>   Use `force = TRUE` to force installation

library(tmap)
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)
#> 
#> Attaching package: 'tmap'
#> The following object is masked from 'package:datasets':
#> 
#>     rivers

tmap_options(check.and.fix = TRUE)

tm_shape(invalid_geometry) +
  tm_fill()
#> Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 9 crosses edge 26

[Nowosad]

Hi @neon-ninja -- see https://www.branchtwigleaf.com/shinyapps/make-valid-geom/. In general -- there is no one way to always fix the geometry issues, and many different geometry issues are possible. Try running sf::sf_use_s2(FALSE) before making the map. If that does not help -- I would suggest you to try to fix the geometry of your data either using the tools shown on the shiny website or manually.

[neon-ninja]

Yes, sf::sf_use_s2(FALSE) does fix it. I suppose if one got the error message "Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 9 crosses edge 26", and put that into their search engine, they might end up on either this page or a page in sf's issue tracker (such as r-spatial/sf#1902) and figure out that they need to set sf::sf_use_s2(FALSE) in addition to tmap_options(check.and.fix = TRUE). This doesn't make for a great developer experience in my opinion though - ideally tmap_options(check.and.fix = TRUE) would work on it's own. It violates the principle of least surprise that tmap_options(check.and.fix = TRUE) doesn't fix invalid geometries.

[Nowosad]

Thanks for your perspective, @neon-ninja.

@mtennekes would it be possible to do something like this in tmap, when tmap_options(check.and.fix = TRUE): 1. check the geometries, 2. if works -- make a map; if fails, try to use st_make_valid(), 3. check the geometries again, 4. if works -- make a map; if fails, switch to the opposite sf_use_s2(), 5. check the geometries again, 6. if works -- make a map; if fails try to use st_make_valid(), 7. check the geometries again, 8. if works -- make a map; if fails return a more descriptive error message?

[mtennekes]

Great idea!

[anjelinejeline]

FYI On 22 January 2024 this is still an issue .. neither tmap_options(check.and.fix = TRUE) or st_make_valid() are useful ...
I solved it by switching to sf::sf_use_s2(FALSE)

[mtennekes]

Implemented using @Nowosad's suggestion. Code review welcome https://github.com/r-tmap/tmap/blob/master/R/check_fix.R

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.3, PROJ 9.2.0; sf_use_s2() is TRUE

temp = tempfile()
tdir = tempdir()
download.file("https://github.com/user-attachments/files/16969923/geoBoundariesCGAZ_ADM2_invalid_polys.rds.zip", temp)
unzip(temp, exdir = tdir)
x = readRDS(file.path(tdir, "geoBoundariesCGAZ_ADM2_invalid_polys.rds"))

sf::sf_use_s2()
#> [1] TRUE

tm_shape(x) +
    tm_polygons() +
    tm_check_fix()
#> The shape object "x" is invalid. Trying to fix it...
#> Shape x has been fixed with s2 = FALSE. If the map doesn't look correct, please run sf::sf_use_s2(FALSE) before running the tmap code again.

sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

tm_shape(x) +
    tm_polygons() +
    tm_check_fix()

^{Created on 2024-09-11 with reprex v2.1.0}

[anjelinejeline]

Implemented using @Nowosad's suggestion. Code review welcome https://github.com/r-tmap/tmap/blob/master/R/check_fix.R

library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.3, PROJ 9.2.0; sf_use_s2() is TRUE

temp = tempfile()
tdir = tempdir()
download.file("https://github.com/user-attachments/files/16969923/geoBoundariesCGAZ_ADM2_invalid_polys.rds.zip", temp)
unzip(temp, exdir = tdir)
x = readRDS(file.path(tdir, "geoBoundariesCGAZ_ADM2_invalid_polys.rds"))

sf::sf_use_s2()
#> [1] TRUE

tm_shape(x) +
    tm_polygons() +
    tm_check_fix()
#> The shape object "x" is invalid. Trying to fix it...
#> Shape x has been fixed with s2 = FALSE. If the map doesn't look correct, please run sf::sf_use_s2(FALSE) before running the tmap code again.

sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

tm_shape(x) +
    tm_polygons() +
    tm_check_fix()

Created on 2024-09-11 with reprex v2.1.0

This works thank you very much @mtennekes

[anjelinejeline]

@mtennekes Hi!!
It has been a long time since we talked about this issue .. I tried to plot another shapefile from the GADM but this time around the sf::sf_use_s2(FALSE) did not work .. haven't you fixed this problem ?

all_countries=readRDS("Countries.RDS")


all_countries_simple=all_countries |> st_simplify(dTolerance = 0.15)

sf_use_s2(FALSE)
tm_shape(all_countries_simple)+
  tm_borders(col='black') +
  tm_check_fix()

[mtennekes]

Which shapefile exactly? And what the output and where there any messages/warnings/errors?

[anjelinejeline]

I just sent it to you by email .. it is too heavy to be uploaded here

[mtennekes]

580 MB is indeed a bit large :-)

Unfortunately this is beyond tmap: also plotting with sf doesn't work well. Probably the spatial object is not well defined. What I did:

library(sf)

x = readRDS("~/Downloads/Countries.RDS")
isv = sf::st_is_valid(x)

table(isv) # 3 were FALSE

ids = which(!isv)

x2 = x[ids,]

rm(x)
gc()

sf::sf_use_s2(FALSE)
plot(x2)
sf::st_is_valid(sf::st_make_valid(x2))

sf::sf_use_s2(TRUE)
plot(x2)
sf::st_is_valid(sf::st_make_valid(x2))

Both draw a map with the horizontal lines, probably due to a 180/-180 longitude crossover.

Also making it valid doesn't fix this artefact.

x3 = sf::st_make_valid(x2)
sf::st_is_valid(x3) # returns all TRUE
plot(x3)

The returned plot is:

Note that sf draws a map per variable, in this case the first 9.

[anjelinejeline]

Hi @mtennekes thank you for looking into it ..
So what should I do to fix this issue ?
Is there anything I could do ?
Thanks
Angela

[nickbearman]

@anjelinejeline The GADM website sometimes has shapefiles with this issue. I'd recommend you try https://geoboundaries.org/ for a shapefile of the same (country/location) area and see if that works.
Probably there is something wrong with the GADM shapefile, perhaps the geometry, and can be tricky and time consuming to diagnose!

[anjelinejeline]

@nickbearman I initially wanted to use geoboundaries but I am not comfortable with its classification .. for instance ADM2 in Italy does not correspond to the right one instead it corresponds to ADM1 as correctly reported in GADM.
I need a world shapefile of adm2 divisions and GADM is the only one I found reflecting the truth .. any other suggestions?

[mtennekes]

Also check out the giscoR package. I've discovered it recently during teaching, and seems pretty good at first sight.

[anjelinejeline]

@mtennekes unfortunately it does not contain adm2 at the global level

[nickbearman]

https://www.naturalearthdata.com/ might also be worth a look