What percentage of polygon x is in raster pixel y, efficiently

I feel this is a competition worthy problem, due to its size.

I would like to weight a set of raster pixels to the attributes of a set of polygons. To do this I need to know the proportion of each polygon in each raster pixel.

I'd also like to keep IDs of both the pixel and the polygon on the final output.

Sometimes the polygons are larger than the pixels, sometimes smaller. Pixels with a value of 0 can be excluded.

The polygons I am working with are Statistical Area Level 1 (SA1) ASGS Ed 2011 and they can be found at the bottom of this page:

http://www.abs.gov.au/AUSSTATS/[email protected]/DetailsPage/1270.0.55.001July%202011?OpenDocument

The raster I am working with is here:

http://www.abs.gov.au/AUSSTATS/[email protected]/DetailsPage/1270.0.55.0072011?OpenDocument

I believe these are GDA94, the polygons are epsg:4283, the raster is epsg:3577

As you'll notice these are both large, which is why efficiency is important.

Any open source solution (R, QGIS etc.) would be best, although if there is an efficient solution in ArcGIS or MapInfo it would be interesting to know of it.

• what are you trying to do here, is it to estimate the population of each polygon from the raster cells it overlaps, by taking a fraction of each pixel's value according to the area of overlap? Aug 5, 2015 at 11:22
• The two projections are different. The raster is in Albers Equal Area (`+proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs`) and the shapefile is in Lat/Long: `+proj=longlat +ellps=GRS80 +no_defs`. Your choice of which projection you map the two datasets together will impact your answer. You need to project one to the other -- which would it be? Aug 7, 2015 at 17:39

from my understanding of your question, and looking at the data, you want to assign a population to each polygon. In pseudocode, something like

``````for each polygon
popn <- 0
for each intersection(polygon,pixel outline)
fraction <- calculate area of intersection (as percentage of pixel outline value)
popn <- popn + (population*fraction)
``````

To do this, you can use QGIS and PostgreSQL/PostGIS

In QGIS,

• convert the raster to polygons. Make sure you check the 'DN' box, so the pixel values are copied over. Note that this won't be a 1:1 mapping between pixels and polygons - contiguous pixels with the same value will form an area. More on this later...

• select all polygons with a value of >0, and save selection as a new shapefile. This will speed things up by removing most of Aus, which is largely empty ;)

Zooming in, you should now see something like this

Note that the irregular shapes have a single value - e.g. 3. That means that each pixel in the shape has a value of 3, so if the shape has the area equivalent to 5 pixels, that shape actually has a population of 15.

• import this shapefile into postgres using shp2pgsql, as table 'pixels'

• do the same for the statistical areas, as table 'areas'

I used epsg:3577, all that matters is that both are the same.

You can now use postgis to get the intersection of each area polygon with each 'pixel area' polygon, and sum up the values for each polygon of

``````(a/b) * (pop*(b/p))
``````
• where a is the area of the intersection with the pixel area
• where b is the area of the pixel-area (which maybe one pixel, or several)
• where p is the area of a single pixel (constant, 1 km^2)
• where pop is value (population) of the pixel-area

that should simplify down to

``````pop * (a/p)
``````

So, something like

``````create table overlap as (
select
a.sa1_main11,
p.DN * (ST_Area(ST_Intersection(a.geom,p.geom))/1000000.0) as weighted_pop,
ST_Intersection(a.geom,p.geom)
from
pixels as p,
areas as a
where
a.geom && p.geom
limit 10000
);
``````

Bringing this back into QGIS, and colourising each intersection at random, you can see how the areas have been parcelled up.

You should be able to see that the populations add up to the same values as the original.

You can now be able to do a query to aggregate the weighted populations on the sa1_main11 field. You can join your original area shapefile to the result of this query.

obviously, take off the limit when you want to do the whole thing. I did this so I could get a result in a couple of minutes to check... but it got in most of Western Australia.

• Thanks Steven Kay - You’ve provided such a good, thorough answer that I think you might have scared off any other contenders! I’m actually weighting the population of each grid square (known) to the sub populations of the polygons (known) – but it’s basically the same as any large scale reweighting. There are some nuances/complexity in the actual weighting that’s not directly relevant to this post.
– K .
Aug 24, 2015 at 1:12

Just in case anyone has come here for a problem like this one here are a few notes:

• epsg: 3577 (as Steven Kay used) is the correct spatial reference in this case (Australia) because it preserves area. Spatial references that don’t exactly preserve area (like epsg:4283) will return a slightly wrong area. Everything will need to be re-projected to the chosen area preserving projection (thanks aaryno for pointing that out.)

• Contiguous pixels with the same value must not form an area - because although the unweighted population is the same the weighted population will be different.

• I think that the correct functions in PostGIS are geography (spherical) rather than geometry (Cartesian).

• The advantage of using geometry is that it is much less complex and thus quicker to run – to the point where it’s feasible. At least under my system constraints, geography is not feasible, even when limited to 10,000, let alone the 500,000 odd required.

• The solution I’ve used is to calculate the overlap based on un-projected geometry as per Steven Kay’s answer. Then as a second step to re-project the resulting polygons into ‘geography’ and calculate the area.