Error in gUnion() function for specific shapefiles

[PieterjanVerhelst]

For the project 2015 PHD VERHELST EEL I tried to combine two shapefiles, one are lines (Vhag.shp) and the other is a polygon (ws_bpns.shp). Note that this is to calculate distances between stations from the Belgian Lifewatch network. I now use other shape files than before, because the 'new' ones have a better fit to the actual rivers and contain more tributaries (= more useful in the future when tracking network potential changes or expands). However, when running gUnion() this gives an error (line 33 - 45):

> study.area <- gUnion(vhag, ws_bpns)
ws_bpns is invalid
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  : 
  TopologyException: Input geom 1 is invalid: Self-intersection at or near point 567066.46833233652 5694796.7879214734 at 567066.46833233652 5694796.7879214734
In addition: Warning messages:
1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point 567066.46833233652 5694796.7879214734
2: In gUnion(vhag, ws_bpns) :
  Invalid objects found; consider using set_RGEOS_CheckValidity(2L)

For other projects I would also like to apply gUnion() on the following shapefile sets:

• Vhag.shp + ws_bpns.shp
• Vhag.shp + ws_bpns.shp + LowCountries_Water_2004.shp
• Vhag.shp + ws_bpns.shp + isabella_kanaal
• ws_bpns.shp + WATERDEEL_VLAK.shp

Note that I use combinations per project, as I think creating a big shapefile covering Belgium and The Netherlands may in some cases be requiring too much calculation time for what I need. Unless you disagree, than we can make a big Belgium-Netherlands map (= Vhag.shp + ws_bpns.shp + LowCountries_Water_2004.shp + WATERDEEL_VLAK.shp)

All shapefiles are uploaded in the folder \data\Belgium_Netherlands under the branch #69 . Hence, I would only address this issue after that branch has been merged with the master branch.

[PieterjanVerhelst]

@damianooldoni just a kind reminder if you have an idea to solve this issue. If this is tackled, I can proceed with the analysis of the data 😄 .

[damianooldoni]

I will fast check whether the solutions discussed in this post work and if yes, which one is the best one.

[damianooldoni]

I find working with sf so chill :smile:

> # transform to sf data.frame
> ws_bpns_sf <- st_as_sf(ws_bpns)
> ws_bpns <- load.shapefile("./data/Belgium_Netherlands/ws_bpns.shp",
+                           "ws_bpns",
+                           coordinate.string)
OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\damiano_oldoni\Documents\INBO\repositories\fish-tracking\scripts\receiver_distance_analysis\data\Belgium_Netherlands\ws_bpns.shp", layer: "ws_bpns"
with 1 features
It has 1 fields
Integer64 fields read as strings:  Id 
> # transform to sf data.frame
> ws_bpns_sf <- st_as_sf(ws_bpns)
> if (isFALSE(st_is_valid(ws_bpns_sf))) {
+   message(glue("Shapefile is invalid. Make it valid first."))
+   ws_bpns_sf <- st_make_valid(ws_bpns_sf)
+   valid_shape <- st_is_valid(ws_bpns_sf)
+   message(glue("Validation succeded: {valid_shape}"))
+ }
Shapefile is invalid. Make it valid first.
Validation succeded: TRUE
> # transform back to SpatialPolygonsDataFrame (sp)
> ws_bpns <- as_Spatial(ws_bpns_sf)
> remove(ws_bpns_sf)

It takes ages to import vhag, I don't know why. But I would apply the same strategy as done above for ws_bpns. If it works, I will put this validation step in a function so that you don't have to write so many lines of code in prepare_dataset.R.

Let me know if it works. Thanks.

[damianooldoni]

vhag seems to have a valid geometry. And so, the commando study.area <- gUnion(vhag, ws_bpns) doesn't throw any error anymore. :+1

@PieterjanVerhelst, please, check this and if you find it working as I found, I will make a branch where I will add the validation as a standard part of the preparation and I will put it in a function. A PR will then follow fast.

[PieterjanVerhelst]

@damianooldoni thanks for the quick reply! Indeed vhag takes a very long time to upload, therefore I created a shapefile with the selection of rivers I need. The file name is pbarn_freshwater.shp and is in the folder \data\Belgium_Netherlands. So if I am correct, you did not made changes to the code (I don't see any commits coming in), but you want me to try the chunk code from above, right?

[PieterjanVerhelst]

update: that code works perfectly 😄 !! Here the result:

Bring on the PR!

[damianooldoni]

Cool! I hope to find a little of time for it tomorrow.