Large shapefile to raster

Question

I have a large shapefile (1 gb) and need to rasterise it. I have tried the following already.

1.) Import it into GRASS using v.in.ogr this failed with the error message: ERROR: G_realloc: unable to allocate 498240036 bytes at break_polygons.c:188

2.) My second idea was to use PostGIS. Import the shapefile, resample it at the x,y location of the grid and then export these points and create a grid from xyz. I successfully imported the shapefiles (polygons and points) but intersecting 1 million polygons with 300k points seems to be very slow. I used the following PostGIS satement, maybe there is room for improvement.

select polygons.land_id,grid.geom from grid,polygons where grid.geom && polygons.geom and within(grid.geom,polygons.geom)

3.) I did also try to use simplify() in PostGIS. But I lost to many small polygons (i.e. some areas that were covered only with small polygons became null).

Any ideas would be greatly appreciated.

Answer 1

You could try gdal_rasterize, although I've not used it with such a large shapefile, so you may have the same issues as you did with GRASS. I reckon something like the following should work (with GDAL >= 1.8.0):

gdal_rasterize -a AN_ATTRIB -l THE_LAYER -a_nodata -9999 -a_srs EPSG:27700 -co TILED=YES -tr 10 10 -ot Float32 src.shp dest.tif

Of course, you will have to play around with some of the options depending on your source shapefile. The most important parameter is -tr which specifies the resolution of a pixel; without it, you may find yourself with a very large raster...

If you want to stick with GRASS, try setting a smaller extent for the rasterization, and split the process up into manageable chunks, then mosaic the rasters into one.

Answer 2

Regarding your first intent, could you try running v.in.ogr command on a machine having more RAM or swap disk space?

If not you can split it in many files an rasterize each of them before merging them again.

Answer 3

If you need more control when doing a rasterization, check out perrygeo's poly_density.py script which uses GDAL under the hood, but can be used for overlapping features or adding conditional evaluations beyond those possible with gdal_rasterize alone.

Answer 4

Like Jack the Ripper, do it by parts. From postgis export slices of the data import in grass and convert.

Answer 5

Here is a guide on how to rasterize a vector table in PostGIS...