How to reproject with gdal_rasterize?

Question

How do I use GDAL tools to rasterise a shapefile, while also reprojecting to a different target coordinate system?

The gdal_rasterize documentation says on-the-fly re-projection is supported in recent versions, but doesn't explain how to perform it. It seems not to document any projection related arguments except -a_srs (which it says is for overriding not reprojecting). Do I need to first create an empty destination file with the target attributes, and if so then what tool do I use to do so?

Answer 1

If the target image already exists then re-projection happens automatically. However, GDAL must recognize the coordinate system of vectors. Check that with ogrinfo. Is source projection is not found (for example because shapefile does not have the .prj part) you can tell it manually by using the -a_srs parameter http://www.gdal.org/gdal_rasterize.html.

Unfortunately it is not possible to create a new target image that is in another coordinate system than vector data. See GDAL issue https://trac.osgeo.org/gdal/ticket/7197.

You can create a blank GeoTIFF to be used as a receiving raster for you vectors with Python as described in https://trac.osgeo.org/gdal/wiki/FAQRaster#HowcanIcreateablankrasterbasedonavectorfilesextentsforusewithgdal_rasterize. I can't discover how to create such an empty raster with any of the GDAL command line tools.

Answer 2

IMO, gdal_rasterize really needs a -t_srs option to reproject vectors. I guess in traditional survey workflows vector data was often delivered in a projection suitable for direct viewing (read: any conformal projection), but these days, at least with more common "consumer-grade" datasets, it's going to be in a geographic CS about 95% of the time. Maps displayed in EPSG:4326 look inherently terrible in most cases, and gdalwarping the raster to solve that problem is usually going to add some artifacts.

Rather than fiddling with extra code to create a blank image it's probably easier to ogr2ogr the vector data into the projection you want, then rasterize the output of that (you also don't have to worry about hand-specifying the extent in this case since gdal_rasterize will be able to do it normally in this case).