Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

"Split Features" and "Cut with Selected Polgyons" turns polygons into linestrings (Geopackage, GeoJSON) #58

Open
aourednik opened this issue Oct 14, 2021 · 4 comments

Comments

@aourednik
Copy link

aourednik commented Oct 14, 2021

"Split Features" and "Cut with Selected Polgyons" functions often turn polygons into linestrings.
Geometry should be strictly enforced to never generate a feature of another kind than the kind of the original (split or cut polygon).
Both GeoPackage and GeoJSON do support multiple geometry types (a mixture of multiplogons and linestrings in a single file and layer), but this does terrible things in other GIS software (when exporting file from QGIS and opening it with something else)

@aourednik
Copy link
Author

aourednik commented Oct 15, 2021

Workaround: for the time being, corrupt Geopackage files (mixing multipolygons and geometry collections of linestrings) can be repaired with the sf package in R. Assuming the path to your geofpackage file is stored in geofile_fullpath

require(sf) # Combines sp, rgeos and rgdal functionality
require(stringr)
require(magrittr)
geodata <- st_read(geofile_fullpath, quiet=TRUE) ;
geodata$geometrytest <- sapply(geodata$geom,function(x){st_dimension(x)})
geodata$geometryempty <- sapply(geodata$geom,function(x){st_is_empty(x)})
geodata <- geodata[!(is.na(geodata$geometrytest) & geodata$geometryempty),]
geodata$geomclass <- sapply(geodata$geom,function(x){class(x)[2]}) 
if (nrow(geodata[geodata$geomclass!="MULTIPOLYGON",]) > 0) {
		for (j in 1:nrow(geodata[geodata$geomclass!="MULTIPOLYGON",])){
			x <- geodata[geodata$geomclass!="MULTIPOLYGON",][j,]$geom
			mypolygons <- st_collection_extract(x) # weirdly, this generates polygons, even though we have linestrings
			geodata[geodata$geomclass!="MULTIPOLYGON",][j,]$geom <- st_multipolygon(mypolygons)%>%st_sfc
		}
}
geodata$geometrytest <- NULL
geodata$geometryempty <- NULL
geodata$geomclass <- NULL
st_write(geodata,geofile_fullpath,append=FALSE)

See also this on StackOverflow: https://stackoverflow.com/questions/57677412/sf-how-to-get-back-to-multipolygon-from-geometrycollection

@bstroebl
Copy link
Owner

bstroebl commented Oct 18, 2021

@aourednik interesting problem, good that you found a workaround. I did only test with Shape and PostGIS where in both cases the data provider enforces the geometry type. What do you suggest how the tool should behave?

  1. ignore such cases and leave the input untouched?
  2. abort process and tell user to fix the input?
  3. run process but tell user afterwards that the result has mixed geometries?

Could you share example data?

@aourednik
Copy link
Author

@bstroebl Thx for reply. An example of a data containing one featur of "geometry collection of linestrings" type (instead of a multipolygon) resulting from cutting is in this commit: https://raw.githubusercontent.com/aourednik/historical-basemaps/044c9fa854d23087fe5f51cd8481ba17a13ccbf8/geojson/world_bc100.geojson

The best solution for my usecase would be to enforce the geometry type in the algorithm itself, i.e., scan the layer for geometry types and use the most frequent one when generating the new, cut, features. This could happen silently or with a warning like "Geometry type of resulting features was determined from the most frequent geometry type in the layer".

I note that your tool does a much better job at keeping geometry types coherent than the preinstalled digitizing tools of QGIS 3. I have the feeling that many QGIS' own and package algorithms remain Shape-oriented, while the software seems to promote the newer GeoPackage format (an excellent step in all other reagrads, since Shape does horrible things to Unicode string fields, among other inconveniences)

@bstroebl
Copy link
Owner

The question is why the cutting results in a linestring. I guess it is because the two polygons are (at least partly) identical. To verify I would need both input layers. If the result of an operation is a linestring one cannot make it a polygon. I can only see the three possible solutions I noted above.
100% agreed regarding the usefulness of shape files.

# for free to join this conversation on GitHub. Already have an account? # to comment
Projects
None yet
Development

No branches or pull requests

2 participants