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

Question

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.

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 1

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.

Use the aggregation tool as follows:

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

Answer 2

using ArcGIS basic you can perform "intersect_analysis" (assuming your grid is made of polygons, otherwise you first need to convert to polygons), then compute the area of the new polygon, then use "summarize table" based on the grid index field (and the land use field as well, optionnally) that will be stored in your new feature class.

the similar workflow can be done in QGIS.

As a remark, with arcgis advanced you can use "tabulate intersection" directly

Answer 3

Thank you both for your very quick answers, which both solve my problem from a process point of view. Using OpenJump is much easier, but requires PC-power. ArcGIS requires licence money and a few more steps before getting there, but less hardware resources.

PS. I dug even further into my long time favorite, GME, and found out that it actually has the functionality I thought it did´nt:

If v is the continuous value specified by the 'field' parameter, c is the area of the (clipped) continuous data polygon falling within the zonal polygon, and A is the total area of the (unclipped) continuous data polygon, then the area weighted mean is calculated as sum(v*c)/sum(c), and the area weighted sum is calculated as sum(v*c/A). In these formulae sum iterates over the set of polygons that overlap that zonal polygon.

Answer 4

The sf package in R has a good solution for this

out <- st_intersection(grid, polygons)

Dr. Grant Williamson at University of Tasmania has a fuller vignette here: https://atriplex.info/blog/index.php/2017/05/24/polygon-intersection-and-summary-with-sf/

A similar question with another approach in R can be found here: https://stackoverflow.com/a/55383754