Friday, 22 December 2017

postgis - Large shapefile to raster


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



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.


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