Wednesday, 21 January 2015

spatial statistics - Summarize the total area of land use polygons within map sheets - how?


I would imagine that at least some of you have encountered similar problems to mine and I´d like to hear how you solved them.


I have a map sheet index (a grid as a shape file) covering a large area. I also have land use polygons covering roughly the same area. Now I would like to calculate how many hectares of land use polygons fall inside each map sheet (i.e. only the portion of the land use polygons that fall within the map sheet polygon should be taken into account when summarizing "land use area"). The picture below might help to see what I mean, map sheets in blue lines, land use parcels in light green.


enter image description here


Earlier the easiest way to do this for me was to use Hawth´s Tools polygon-in-polygon analysis (http://www.spatialecology.com/htools/polypolyanalysis.php), it did exactly what I wanted. Hawth´s Tools is now replaced by GME and the polygon-in-polygon analysis tool in GME no longer functions the same way as the old one. The new tool only sums up the total area of all land use parcels that intersect each map sheet and writes that as the result. No good for me.


I could use the map sheet lines to cut the land parcel dataset into small bits and then re-calculate the area for all polygons, but that´s not that tempting: I have 26 000 map sheets and 1,2 million land use parcels.


I´ve tried this in ArcGis Basic but no luck there. OpenJump, Qgis or gvSIG do not seem to have a solution either. Do anyone of you?




Answer



OpenJUMP does have a solution. You need the Plus version which includes the Aggregation plugin (Plugins - Analysis - Aggregation).


EDIT I took execution time from a test where I used similarly sized datasets. I created a polygon layer with 1.2 million polygons and a polygon grid to present 26000 map sheet rectangles. Computing the parcel area per mapsheet took 55 seconds with a good desktop. More than 5 GB of memory is needed as well as 64 bit java.


Following images show how this tool is used. Layer "New" is the polygon layer, "New (2)" represents the map sheet layer. The area of the polygon is 39123 units.


enter image description here


Use the aggregation tool as follows: enter image description here


Sum of intersected areas is moved as attribute into the map sheet layer.


enter image description here


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