Thursday, 17 September 2015

qgis - Getting error when clipping shapefile to shapefile using ogr2ogr?


I want to clip a shapefile to a shapefile. The input shapefile is a 1 m vector grid generated in QGIS. The tile I'm clipping to is just a simple polygon (the geometry is contained within the 1 m vector grid).


I put into the command line:



ogr2ogr -overwrite -t_srs EPSG:27700 -clipsrc ./Crewe_boundary_tile_grid_ID__0.shp ./Crewe_boundary_1m_grid_clipped_ID__0.shp ./Crewe_boundary_1m_grid.shp



However I get the error message:




ERROR 1: Attempt to write non-polygon (LINESTRING) geometry to POLYGON type shapefile. ERROR 1: Unable to write feature 500 from layer Crewe_boundary_1m_grid.


ERROR 1: Terminating translation prematurely after failed translation of layer Crewe_boundary_1m_grid (use -skipfailures to skip errors)



Yet if I try and run this for a different shapefile, i.e. a simple rectangular polygon, it works. I have also used exactly the same command on a 1000 m grid and that works. Also, if I do 'ogrinfo' then each of the input files reads as a polygon.


How do I fix this?



Answer



The ogr2ogr clip operation creates polygons but for some reason also linestrings to to South and East edges of the area. Pink lines below show those 100 linestings (one is selected). Because of mixed geometrytypes the result cannot be saved into shapefile and therefore the error. I am not sure if this is an intended behaviour of ogr2ogr clipping or a bug. Anyway in your case it is nothing to be worried.


enter image description here


You can safely select to write a polygon shapefile and skip the errors. Using explicit -nlt is a good habbit or otherwise the new shapefile will take the geometrytype of the first feature that is met and that is not necessarily the one that is wanted.




ogr2ogr -f "ESRI shapefile" -clipsrc clip.shp output.shp input.shp -nlt POLYGON -skipfailures



Next time if you want to see and understand what happens you can make conversion into GML format which supports mixed geometries. Better way for my mind is to convert into Spatialite with a generic geometry type and do the rest of the research with SQL. The first steps could be



ogr2ogr -f "SQlite" -dsco spatialite=yes -clipsrc clip.shp output.sqlite input.shp -nlt GEOMETRY



create table lines as
select * from crewe_10m_grid
where geometrytype(geometry)='LINESTRING';


create table polygons as
select * from crewe_10m_grid
where geometrytype(geometry)='POLYGON';

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...