Tuesday, 19 May 2015

rasterization - Area-weighted rasterising of polygon data


I'm looking for any GIS tool to do detailed rasterising of polygons. I've checked GDAL's gdal_rasterize command (and API), and ArcGIS's Polygon to Raster tool, but they only rasterise 1 attribute per cell, not an area-weighted aggregated attribute. Let me describe with an example. Consider that I need to rasterise values from polygons into a 100 m resolution raster. Here are three polygons contained in one pixel:


polygons in one pixel


Where the areas for each class are 4000, 2500, and 3500 m², or a weighting of 0.4, 0.25 and 0.35. If the values of the attribute that I need to rasterise are 12, 42 and 8, then their weighted values are 4.8, 10.5 and 2.8. The sum of the weights, 18.1, is the value that I would like to rasterise, since it is representative for the whole pixel. Most other rasterisation routines I'm aware of will choose 12, since that is the value of one of the polygons that is either at the centre of the cell or has the maximum area. 12 is not acceptable for my analysis, since it will skew the mass-balance for the whole raster, which is unacceptable.


I'm looking to see if there are better tools for rasterising. My current approach is to perform a vector analysis at the pixel level, using the intersection of polygons and each pixel, then combine the results into a raster. This is a slow process, as I have 80275 rows (each pixel).



Answer



The approach mentioned by @whuber is the quickest, even if it isn't 100% perfect. Here are the steps using GDAL tools (at least version 1.10) from the command line.



The finer resolution raster needs to be fine enough to reasonably capture the detail of the vector data, which depends on the vector source. The fine resolution should be a multiple of the coarse resolution, and the extents should also be a multiple of the coarse raster. See more about gdal_rasterize.



E.g. rasterise myattr from polygons.shp into a new 10 m resolution raster vals_10m.tif with specified extents defined by the -te window, which are 21000 m by 29500 m, so the resulting raster should be 2100x2950:


gdal_rasterize -te 888200 749000 909200 778500 -tr 10 10 -a vals -ot Float32 \
-a_nodata -999 -init -999 polygons.shp vals_10m.tif

first



This can be done using the gdalwarp tool using the average resampling method to computes the average of all non-NODATA contributing pixels. The extents for both fine and coarse resolution rasters should be the same, if the correct extents and resolutions are chosen in the first step.


E.g. resample the fine resolution (10 m) vals_10m.tif to a coarse resolution (100 m) vals_100m.tif that is 210x295:


gdalwarp -tr 100 100 -r average -dstnodata -999 vals_10m.tif vals_100m.tif


second


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