2

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?

  • Have you tried FME? safe.com Feature Manipulation Engine can tranform your data with a buffer and union in a quick translation – Mapperz Aug 16 '16 at 15:38
  • I have only tried what is stated above. I do not know what FME is... – J.A.Cado Aug 16 '16 at 15:39
  • Updated the original comment with details. – Mapperz Aug 16 '16 at 15:44
  • Even if you drop your PostgreSQL connection the query should still be running until you cancel it, also your query is buffering 0.05 degrees (edit. saw your point of it being about 5.5km, you should read: en.wikipedia.org/wiki/Decimal_degrees) – HeikkiVesanto Aug 16 '16 at 15:56
  • Have you tried the 64 bit version of QGIS? – klewis Aug 16 '16 at 17:14
3

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;

  • What EPSG would your recommend for a global Cartesian coordinate system? Thus far for anything global I've been using an equal-area proj4 string with +proj=cea but I'm wondering if you can suggest something better. – srha Oct 31 '17 at 22:24
  • 1
    Actually for storing world data I would recommend a geographic coordinate system (like EPSG 4326). However when it comes to projecting it, I often use EPSG 3857, if precision is not a matter (this is the projection on the sphere!). It is good for overlaying your data on popular basemaps (Google, ESRI and so on), so you don't need to reproject on the fly. Another variant is EPSG 3395 (WGS 84 / World Mercator), which is projection on the world ellipsoid, and it gives more precise measurements. – Yaroslav Nov 1 '17 at 15:01
2

I've just tried to use the QGIS Vector Geometry tools-> Fixed distance buffer on GSHHS_I_L1 and it completed in less than a second.

So I suspect you are using a much higher resolution version of the shapefile (There is no need if you are studying the whole world) or you are on a much smaller machine than I am.

  • 1
    Although you have a very valid point, and indeed I did pick the largest shape file when there is absolutely no need to, I did not accept your answer because it doesn't answer my question. But again....it's very very good point you make. So thanks :D – J.A.Cado Aug 22 '16 at 15:21
2

I doubt you can get much faster than PostGIS. You can script it to run in batches (in your favourite shell or language), so the number of total features won't be a problem. Just create the table and update it in steps, only handling a few features at a time (using LIMIT and OFFSET).

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.