neon-ninja commented 2 years ago

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 commented 2 years ago

Have you tried to set st_use_s2(FALSE)?

neon-ninja commented 2 years ago

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 commented 2 years ago

Sorry - I meant sf_use_s2(FALSE)

neon-ninja commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

Looks like someone already raised this issue over in the sf repo - https://github.com/r-spatial/sf/issues/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 commented 2 years ago

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 commented 2 years ago

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 https://github.com/r-spatial/sf/issues/1780#issuecomment-911729141)?

eblondel commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

@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 commented 1 year ago

[https://res-h3.public.cdn.office.net/assets/mail/file-icon/png/zip_16x16.png]Tmap data files.RData.ziphttps://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)

@pennybeaverhttps://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 GitHubhttps://github.com/r-tmap/tmap/issues/606#issuecomment-1706329872, or unsubscribehttps://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 commented 1 year ago

Tmap data files.RData.zip

zip file attached :)

Nowosad commented 1 year ago

@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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 https://github.com/r-spatial/sf/issues/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 commented 1 year ago

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 commented 1 year ago

Great idea!

anjelinejeline commented 8 months ago

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 commented 1 week ago

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 commented 1 week ago

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

r-tmap / tmap

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

Breeding populations locations

obtain bathymetry data for map