Rasterizing polygons with the function "gdal_rasterize" in R

Question

I'm using the function gdal_rasterize in R (package gdalUtils) to convert a polygon shapefile into a raster. Here is my code:

r.test <- gdal_rasterize(src_datasource="df.shp", dst_filename="test.tif", a="class_raster", output_Raster=TRUE))

"class_raster" are numeric values which should be transferred in the raster.

But I have this error message:

Error in .local(.Object, ...) : 
  `C:\Users\Documents\test.tif' does not exist in the file system,
and is not recognised as a supported dataset name.


In addition: Warning message:
running command '"C:\OSGeo4W64\bin\gdal_rasterize.exe" -a "class_raster" "df.shp" "test.tif"' had status 1 
Error in .rasterObjectFromFile(x, objecttype = "RasterBrick", ...) : 
  Cannot create a RasterLayer object from this file. (file does not exist)

Here is a description of the shapefile:

> df
class       : SpatialPolygonsDataFrame 
features    : 82077 
extent      : 265380.5, 462477.9, 4984446, 5081957  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 1
names       : class_raster 
min values  :            1 
max values  :            8

Update:

I created a file .tif before to call the function gdal_rasterize as follows:

  r.test <- raster(extent(df), resolution = 30) 
  projection(r.test) <- CRS("+proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

> r.test
class       : RasterLayer 
dimensions  : 3250, 6570, 21352500  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 265380.5, 462480.5, 4984457, 5081957  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +towgs84=0,0,0 


> test <- writeRaster(r.test, filename="test.tif", format="GTiff", overwrite=TRUE)
Warning message:
In .local(x, filename, ...) : all cell values are NA

> gdal_rasterize(src_datasource="df.shp", dst_filename="test.tif", a="class_raster", output_Raster=TRUE)
class       : RasterBrick 
dimensions  : 3250, 6570, 21352500, 1  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 265380.5, 462480.5, 4984457, 5081957  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-73.5 +k=0.9999 +x_0=304800 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
data source : H:\test.tif 
names       : test 

But I still have the message:

Warning message:
running command '"C:\OSGeo4W64\bin\gdal_rasterize.exe" -a "class_raster" "df" "test.tif"' had status 1 

Here is a description of the raster with the function gdalinfo

> gdalinfo("H:/test.tif")
 [1] "Driver: GTiff/GeoTIFF"                                                     "Files: H:/test.tif"               
 [3] "Size is 6570, 3250"                                                        "Coordinate System is:"                                                    
 [5] "PROJCS[\"unnamed\","                                                       "    GEOGCS[\"NAD83\","                                                    
 [7] "        DATUM[\"North_American_Datum_1983\","                              "            SPHEROID[\"GRS 1980\",6378137,298.2572221010002,"             
 [9] "                AUTHORITY[\"EPSG\",\"7019\"]],"                            "            TOWGS84[0,0,0,0,0,0,0],"                                      
[11] "            AUTHORITY[\"EPSG\",\"6269\"]],"                                "        PRIMEM[\"Greenwich\",0],"                                         
[13] "        UNIT[\"degree\",0.0174532925199433],"                              "        AUTHORITY[\"EPSG\",\"4269\"]],"                                   
[15] "    PROJECTION[\"Transverse_Mercator\"],"                                  "    PARAMETER[\"latitude_of_origin\",0],"                                 
[17] "    PARAMETER[\"central_meridian\",-73.5],"                                "    PARAMETER[\"scale_factor\",0.9999],"                                  
[19] "    PARAMETER[\"false_easting\",304800],"                                  "    PARAMETER[\"false_northing\",0],"                                     
[21] "    UNIT[\"metre\",1,"                                                     "        AUTHORITY[\"EPSG\",\"9001\"]]]"                                   
[23] "Origin = (265380.523957236320000,5081957.489292523800000)"                 "Pixel Size = (30.000000000000000,-30.000000000000000)"                    
[25] "Metadata:"                                                                 "  AREA_OR_POINT=Area"                                                     
[27] "Image Structure Metadata:"                                                 "  COMPRESSION=LZW"                                                        
[29] "  INTERLEAVE=BAND"                                                         "Corner Coordinates:"                                                      
[31] "Upper Left  (  265380.524, 5081957.489) ( 74d 0'28.08\"W, 45d52'34.80\"N)" "Lower Left  (  265380.524, 4984457.489) ( 73d59'59.97\"W, 44d59'56.44\"N)"
[33] "Upper Right (  462480.524, 5081957.489) ( 71d28' 9.73\"W, 45d51'33.89\"N)" "Lower Right (  462480.524, 4984457.489) ( 71d30' 2.06\"W, 44d58'57.36\"N)"
[35] "Center      (  363930.524, 5033207.489) ( 72d44'39.21\"W, 45d26'10.68\"N)" "Band 1 Block=6570x1 Type=Float64, ColorInterp=Gray"                       
[37] "  Min=1.#IO Max=-1.#IO "                                                   "  Minimum=1.#IO, Maximum=-1.#IO, Mean=1.000, StdDev=1.000"                
[39] "  NoData Value=-1.6999999999999999e+308"                                   "  Metadata:"                                                              
[41] "    STATISTICS_MAXIMUM=-1.#INF"                                            "    STATISTICS_MEAN=1.#SNAN"                                              
[43] "    STATISTICS_MINIMUM=1.#INF"                                             "    STATISTICS_STDDEV=1.#SNAN" 

I have the version 1.0-7 of rgdal.

Answer 1

I cannot recover your error but it seems that the parameters b=... l=... for band and layer are not present and I don't know what is in your data field a=class_raster, may be it is not a numeric type that fits into the raster band, that you want to fill with these data or NA could also be a problem.

Here a little how to, I use to test the gdal_raster in GNU-R (also new for me ..but interesting), may be you can adopt something to get your stuff running. Normally I use the gdal utility at command line under LINUX.

I've two sources in the test case, a 3d shaded DEM called dem.v2.3d.tif and a water mask called wrgn.shp

I made a screen shot of qgis with a 50% reduced screen size for upload. I zoomed a little bit in, so that you can see the underlaying bathymetry, which should get dark blue in the next step . In the R script I use a copy of the dem.v2.3d.tif called test.tif and will blank the water area with the param set ...b=c(1,2,3), burn=c(0,0,128)... (for a dark blue color).

Here is the script test-rgdal.R:

# load the packages
require("rgdal")
require("gdalUtils")

# setup the datat directory
setwd("~/dev.d/gis-exchange.d/data.d/kjb.d");

# copy an existing 3D shaded DEM to the test TIF
file.copy("dem.v2.3d.tif", "test.tif", overwrite=TRUE);


# plain GDAL command to make the water area dark BLUE
# gdal_rasterize -b 1 -b 2 -b 3 \
#   -burn 0 -burn 0 -burn 128 \
#   -l wrgn wrgn.shp test.tif

# the same procedure in R gdalUtil using gdal_rasterize
r.test <- gdal_rasterize(b=c(1,2,3),
                         burn=c(0,0,128),
                         l="wrgn",
                         src_datasource="wrgn.shp",
                         dst_filename="test.tif",
                         output_Raster=TRUE)
# ask whats in
print(r.test)

#EOF

I get this result as expected..

and this output from the script:

> source('~/dev.d/gis-exchange.d/test-gdal.R')
class       : RasterBrick 
dimensions  : 3120, 2680, 8361600, 3  (nrow, ncol, ncell, nlayers)
resolution  : 9.99811, 9.997878  (x, y)
extent      : 385864.4, 412659.4, 6024035, 6055228  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /home/iaw/dev.d/gis-exchange.d/data.d/kjb.d/test.tif 
names       : test.1, test.2, test.3 
min values  :      0,      0,      0 
max values  :    255,    255,    255 

So every thing seems to be fine when I work under a LINUX amd64 environment, DEBIAN Jessie, R Version 3.1.1, gdalUtils 2.0.7.1. on an T410.

Conclusion: You could make some tests with parameters b=... l=... for band and layer and a burn=... param and later go on to use the a=... param.

Answer 2

I'm not sure which version of gdal you're using, but if it's older than 1.8, you have to create the tiff with a proper extension and resolution before you use it.

http://www.gdal.org/gdal_rasterize.html

Look at:

dst_filename:

    The GDAL supported output file. 
    Must support update mode access. 
    Before GDAL 1.8.0, gdal_rasterize could not create new output files.