3

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?

2

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

  • I am using the PLUS-version but I have´nt dug deep enough to note this particular functionality. It does the trick, I did a test on part of my dataset. But when I try to run the tool on my full dataset, OJ runs out of memory when trying to load all the 1,2 million land use parcels. Thanks for the tip, though, I´ll be sure to use it on smaller datasets. – GISinHelsinki Sep 22 '14 at 12:34
  • 32 bit jre can use 2 GB of memory at maximum. 64 bit jre can use as much RAM as the computer has free but if it is not enough then nothing helps because OpenJUMP runs totally in the memory. You may get the data to fit in memory by deleting all the attributes because for summing up areas you do not need them. If that does not help you must split the original dataset. – user30184 Sep 22 '14 at 12:43
  • Thanks for your info on execution time. It seems my problem really is a hardware issue: I need an updated & more powerful desktop. – GISinHelsinki Sep 22 '14 at 12:49
  • OpenJUMP does not generally use much hardware resources: CPU, disk space, not even RAM. However, because all the vector data must fit into memory the amount of memory reserved for java can become a limiting factor. If java must start to swap memory into disk everything comes very slow or OpenJUMP jams totally. – user30184 Sep 22 '14 at 13:16
  • I still run 32-bit Java on a laptop with 4 GB RAM, so I definitely need to update my workstation. But its good to know that its not all hardware, but rather how you tweak your software to use it that does the trick. And I´m really happy to have found this feature in OJ, its one of my favorite desktop GIS. – GISinHelsinki Sep 22 '14 at 13:23
1

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

  • This does the trick, sort of. Its much more laborious and complicated than doing it with OJ, but it seems I can actually get all 26 000 map sheets processed this way. – GISinHelsinki Sep 22 '14 at 12:51
0

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.

0

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

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.