Friday, 27 February 2015

sql - PostGIS Intersection and Summarise Attributes


I am usually an ArcGIS desktop user, but keen to start using PostGIS more and I have a really big bit of processing to do. Not sure what functions to use, hopefully someone can help.


I have a polygon dataset (several million features), based on a type of landuse/ landcover classfication (20 categories). I have a number of regions in another dataset.


For each of the regions, I would like to know the area of each landcover classfication.


In ArcGIS (if it was a smaller dataset) I would imagine first adding the region to each of the polygons in the attribute table using a join. Then using "summarize" on the table by region and by landcover classification.


Not sure where to start doing this in PostGIS / SQL.


Update:


Wow thanks that has been a huge help.


It has been running a long time (44 hours!) and now I get:



NOTICE:  TopologyException: found non-noded intersection between LINESTRING (coords
edited) and LINESTRING (coords edited) at (position edited)
ERROR: GEOS Intersection() threw an error!

********** Error **********

ERROR: GEOS Intersection() threw an error!
SQL state: XX000

I assume this is a problem in the original data - just a case of reviewing the original data or can I first check the topology some how for the whole data? Is there something about accepting certain errors / processing tolerances?




Answer



Assuming you have the following table layout


landcover(id,type,the_geom)
region(id,name,the_geom)

Area values per landcover type per region can be calculated using


SELECT r.id, r.name, r.the_geom, l.type, 
Sum(ST_Area(ST_Intersection(l.the_geom,r.the_geom)))
FROM landcover l,
region r

WHERE l.the_geom && r.the_geom
GROUP by r.id, r.name, r.the_geom, l.type

ST_Intersection is used to account for landcover polygons that are only partially within a region.


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