[GIS] How to get a bbox selection of huge shapefile as GeoJSON

fionageojsonpostgisshapefilevector

I have a 4.29 GB Shapefile from which I want to extract features as GeoJSON, defined by a bounding box.

This should be as fast as possible, because I want to use it in a web application.

My first attempt at this was to convert to SQL using shp2psql and import it into PostGIS. I created a spatial index (gist) and did a full ANALYZE VACUUM on the table. The VM has 4 GB of memory.

The following query works, but took 623 seconds(!):

SELECT
    ST_AsGeoJSON(geom)AS geometry,
FROM
    shapes
WHERE
    geom && ST_GeogFromText(
        'POLYGON((5.117705762386322 52.09243663877818,5.118895322084427 52.09243663877818,5.118895322084427 52.09314276719956,5.117705762386322 52.09314276719956,5.117705762386322 52.09243663877818))'
    )
AND ST_Intersects(
    geom,
    ST_GeogFromText(
        'POLYGON((5.117705762386322 52.09243663877818,5.118895322084427 52.09243663877818,5.118895322084427 52.09314276719956,5.117705762386322 52.09314276719956,5.117705762386322 52.09243663877818))'
    )
)

An equivalent using ST_DWithin took even longer.

I decided to try an Amazon RDS PostgreSQL 9.3.1 instance with 68 GB of memory and 26 CPU's. Query time dropped to around 30 seconds. Faster but still not fast enough.

Preprocessing will take forever, even on Amazon RDS. And it will cost a lot of money.

I also tried Fiona and 'ogr2ogr -spat', but they're about as slow as PostGIS.

What other options exist for quickly getting GeoJSON out of a large Shapefile for a certain spatial extent?

Thanks!!

Best Answer

Since you will go for it, I would like to suggest you some links that may give you some guidance to do it so, here they are:

Definitely you need to adapt them to your needs, but I think it is good to have these links at hand for more people having the same issue you have.

Hope this helps,

Related Solutions

[GIS] Loading geojson via ajax in leaflet using bbox

I hope I understood what You mean. I've changed tho code like this:

var n2k_dh_geojson = L.geoJson(null);

map.on('dragend', function onDragEnd(){
n2k_dh_geojson.clearLayers();

$.ajax({
    type: "GET",
     url: "./mapdata/n2k_dh_geojson.php?bbox=" + map.getBounds().getWest() + "," + map.getBounds().getSouth() + "," + map.getBounds().getEast() + "," + map.getBounds().getNorth(),
    dataType: 'json',
    success: function (response) {
        n2k_dh_geojson = L.geoJson(response, {
            onEachFeature: function (feature, layer) {
            layer.bindPopup(feature.properties.sitename);
            }
        });
        n2k_dh_geojson.addTo(map);
    }
});
});

Now the 'old' data is erased from the map before the 'new' data is drawn, which produces always a little decay in time, which isn't nice or is there a way to cache the old data?

[GIS] Compare only valid polygon geometries in PostGIS

This should skip over invalid geometries.

SELECT
   DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
(CASE WHEN ST_INTERSECTS(t1.geom, t2.geom) IS false THEN '9999.99'
      ELSE ST_HausdorffDistance(t1.geom, t2.geom)
END) as hausdorff
 FROM 
    (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345) as t1,
   (SELECT id, gid, geom FROM old_lyr WHERE ST_IsValid(geom) ) as t2
 WHERE
      ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

You might want to compare teh performance with putting the ST_IsValid check in the outter where. I think both should work but the planner isn't guaranteed to run them in order.

SELECT
   DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
(CASE WHEN ST_INTERSECTS(t1.geom, t2.geom) IS false THEN '9999.99'
      ELSE ST_HausdorffDistance(t1.geom, t2.geom)
END) as hausdorff
 FROM 
    (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345) as t1,
   (SELECT id, gid, geom FROM old_lyr ) as t2
 WHERE
      ST_IsValid(t2.geom)  AND ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

Update to original answer to force the planner's hand.

Here I use ST_Expand && which is just a bounding box check so should be safe with invalid geometries as well. This will materialized the t2 table for those where bounding box is close enough to t2. Using && should be safe with even invalid geometries since invalid geometries should have a bounding box still. This may on downside produce a large worktable with false-positives.

Note the next part I use _ST_DWithin instead of ST_DWithin. The reason is that ST_DWithin is just syntactic sugar for ST_Expand(geom,2) && geom and _ST_DWithin, so no need to do the && again if you've already done it.

 WITH t1 AS (SELECT id, gid, geom FROM new_lyr WHERE gid = 12345),
      t2 AS (SELECT o.id, o.gid, o.geom FROM old_lyr AS o 
       INNER JOIN t1 ON  ST_Expand(t1.geom, 2) && o.geom
             WHERE ST_IsValid(o.geom) )
 SELECT
   DISTINCT ON (t1.gid) t1.gid, t1.id, t2.gid, t2.id,
(CASE WHEN ST_INTERSECTS(t1.geom, t2.geom) IS false THEN '9999.99'
      ELSE ST_HausdorffDistance(t1.geom, t2.geom)
END) as hausdorff
 FROM 
    t1,
    t2
 WHERE
      _ST_DWithin(t1.geom, t2.geom, 2)
ORDER BY
    t1.gid, hausdorff;

Best Answer

Related Solutions

[GIS] Loading geojson via ajax in leaflet using bbox

[GIS] Compare only valid polygon geometries in PostGIS

Related Question