Friday 22 March 2019

qgis - Buffer and union of large shapefile?


I want to create a buffer of 5km around the coastline, and save that as a separate shapefile.


I have a shapefile with almost 180,000 features. It's a line shapefile, representing the coastlines of the world (I got the data here). The shapefile is in CRS (EPSG:4326, WG84).


So far I have been testing three approaches and all are either very slow or crash:


1) QGIS 2.14.2 Essen (with functioning GRASS functions) I go to 'Vector' --> 'Geoprocessing Tools' --> 'Buffer(s)...' I select the shapefile, For buffer distance I use: 0.05 (EPSG:4326 which is about 5.5km) I check the 'Dissolve buffer results' options. This option doesn't seem to work, as QGIS crashes. I expect it has to do with a shortage of memory available. "Error: minidump written to (filepath)" [Screen shot of error message in QGIS[2]


2) PgAdminIII 1.20.0 with the following state I try to create a new table and fill it with the buffer:


CREATE TABLE public.lines_buffer AS SELECT ST_Buffer(geom,0.05) FROM public.lines


But the connection keeps on getting lost and thus breaking the operation. The furthest I got was 2 hours of query time, but that wasn't enough.


3) Python 2.7 (shapely, fiona, descartes) My third attempt, is still in progress. I am trying to loop through the large shapefile and make a buffer. I'm not sure if my code is very efficient. I still have to "dissolve" after the buffer, which would also take a while. I plan on testing shapely.ops.cascaded_union. EDIT: Even running the first 2 items takes very long...[Spyder IDE: 100% memory usage] My current loop just limits the first 3 items in the fiona.collection Credits go to Tom Macwright


def readshape(pathfolder, filename):
filepath = os.path.join(pathfolder, filename)
with collection(filepath) as Input:
schema = Input.schema
with collection('buffer_shape.shp', 'w', 'ESRI Shapefile', schema) as output:
for i in range(0, 3):
output.write({

'properties': {
'source': Input[i]['properties']['source']
},
'geometry': mapping(shape(Input[i]['geometry']).buffer(0.05))
})

Does anyone have any tips, which software is the quickest and most stable option?



Answer



PostGIS should work for your task but:




  1. Reproject your data into some Cartesian coordinate system (like Web Mercator) because computations on the WGS84 ellipsoid (geographic coordinates) in PostGIS take longer time. Also for geography ST_Buffer may not behave as expected if object is sufficiently large that it falls between two UTM zones or crosses the dateline, which is true in your case (PostGIS docs).

  2. Create a spatial index for your features:



CREATE INDEX IF NOT EXISTS lines_idx ON public.lines USING GIST (geom);




  1. Run your query with buffer distance in meters (if you chose a metric CRS):




CREATE TABLE public.lines_buffer AS SELECT ST_Buffer(geom,5500) FROM public.lines;



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...