3

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.

5

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

screen shot qgis

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

resulting image from the R-Script

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.

  • Thank you very much huchfinn for your help. I don't know how to specify colors for the different land cover types in my shapefile "df". In "df", I have 8 land cover types that are classified from 1 to 8 (column "class_raster" with numeric values). Here, "1" corresponds to forests and should be drawn in green (burn=c(0,255,0)), "2" are agricultural fields (burn=c(255,255,0) for a yellow color), etc. Should I build a raster with rgdal_rasterize for each land cover type with c(1,2,3) and burn=c(0,255,0) for forests or c(1,2,3) and burn=c(255,255,0) for fields? – Pierre Feb 11 '16 at 17:51
  • Color isn't nessesary if you use single attributed approach, I guess, if your data field is an integer or a byte you have to design a single band int32 or byte geotif corresponding to the field type. You can use the shplib or ogrinfo to investigate the nature of field types. I recommend you to use the gdal command line tools to investigate your stuff and later adopt the workflow to R. The colouring could later be done by the raster symbolizer. – huckfinn Feb 11 '16 at 17:54
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.
  • Thank you very much huckfinn for your answer. I created the tiff before to call the function gdal_rasterize as follows: test <- writeRaster(r.test, filename="test.tif", format="GTiff", overwrite=TRUE). All raster cells are NA. There is no more the error message: C:\Users\Documents\test.tif' does not exist in the file system, and is not recognised as a supported dataset name but I still have the message: running command '"C:\OSGeo4W64\bin\gdal_rasterize.exe" -a "class_raster" "df.shp" "test.tif"' had status 1 . – Pierre Feb 10 '16 at 23:30
  • Do you build a normal tiff or a geotiff? – huckfinn Feb 10 '16 at 23:33
  • Geotiff. I followed the example in this link: http://www.inside-r.org/packages/cran/raster/docs/writeRaster – Pierre Feb 10 '16 at 23:49
  • Can you check, that the TIFF you've created has the extention like the stuff in your shape file ..for example with gdal utils like gdalinfo? Can you post your code? Mention ..with a proper extention and resolution before you use it. – huckfinn Feb 11 '16 at 0:01
  • Thank you very much huckfinn for your answer. I updated my question by adding my code for the creation of the tiff. Let me know if any precisions need to be made (I am a beginner in GIS and R). – Pierre Feb 11 '16 at 15:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.