1

I have a regular grid of polygons and want to clip them using ogr2ogr. While using sqlite files as destination I stumbled over some wierd behavior:

When clipping polygons exactly at their borders, a mix of geometry types is created so that I end up with points, lines and polygons while I only would expect polygons. When using an ESRI Shapefile, I get errors, thus I add -skipfailures. However, eventually, the Shapefile will only contain one polygon. The figure at the bottom shows the example data.

So my Question is: Is there any possibility to define or select the output geometry type from the clipping? I stumbled over a similar post. Adding to this, would there be any solution, how to do this in one ogr2ogr statement?

I prepared a reproducible example:

#! /bin/bash
echo "ID, GM" > poly.csv
echo '1, "POLYGON((4249000 2863000, 4250000 2863000, 4250000 2862000, 4249000 2862000, 4249000 2863000))"' >> poly.csv
echo '2, "POLYGON((4250000 2863000, 4251000 2863000, 4251000 2862000, 4250000 2862000, 4250000 2863000))"' >> poly.csv
echo '3, "POLYGON((4250000 2862000, 4251000 2862000, 4251000 2861000, 4250000 2861000, 4250000 2862000))"' >> poly.csv
echo '4, "POLYGON((4249000 2862000, 4250000 2862000, 4250000 2861000, 4249000 2861000, 4249000 2862000))"' >> poly.csv

# vrt in order for reproducibility
echo '<OGRVRTDataSource>
    <OGRVRTLayer name="poly">
       <SrcDataSource>poly.csv</SrcDataSource>
      <GeometryType>wkbPolygon</GeometryType>
 <LayerSRS>PROJCS["ETRS89 / ETRS-LAEA",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["latitude_of_center",52],PARAMETER["longitude_of_center",10],PARAMETER["false_easting",4321000],PARAMETER["false_northing",3210000],AUTHORITY["EPSG","3035"],AXIS["X",EAST],AXIS["Y",NORTH]]</LayerSRS>
 <GeometryField encoding="WKT" field="GM"> </GeometryField>
     </OGRVRTLayer>
</OGRVRTDataSource>' > poly.vrt

# convert vrt to sqlite
ogr2ogr -f "SQLite" poly.sqlite poly.vrt

# clip twice
ogr2ogr -f "SQLite" -clipsrc 4250000 2860000 4252000 2862000 out.sqlite poly.sqlite
ogr2ogr -f "ESRI Shapefile" -clipsrc 4250000 2860000 4252000 2862000 out.shp poly.sqlite -skipfailures
# check the result
ogrinfo out.sqlite -al 
ogrinfo out.shp -al 

Output:

# the sqlite
....
OGRFeature(poly):1
  id (String) = 1
  POINT (4250000 2862000)

OGRFeature(poly):2
  id (String) = 2
  LINESTRING (4251000 2862000,4250000 2862000)

OGRFeature(poly):3
  id (String) = 3
  POLYGON ((4250000 2862000,4251000 2862000,4251000 2861000,4250000 2861000,4250000 2862000))

OGRFeature(poly):4
  id (String) = 4
  LINESTRING (4250000 2862000,4250000 2861000)
....

# the shp
....
  POLYGON ((4250000 2862000,4251000 2862000,4251000 2861000,4250000 2861000,4250000 2862000))
....

Figure: visualization of the problem.
blue dashed line: exemplary clipping extent
red: location where the point is created
green: location where the lines are created

1

Ok, I found the solution myself: Just add the following flags-nlt POLYGON -dsco SPATIALITE=YES -skipfailures. Obviously the dataset creation option dsco SPATIALITE=YES prevents the driver to insert differing geometry types into the column. Hence, the desired / expected output is met. skipfailures must be added since, the driver will exit on an error.

ogr2ogr -f SQLite -nlt POLYGON -dsco SPATIALITE=YES -skipfailures -clipsrc 4250000 2860000 4252000 2862000 out.sqlite poly.sqlite 

Unfortunately, I did not yet find any source that would explicitly explain this specific behavior.

Your Answer

By clicking "Post Your Answer", you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.