0

I am trying to rasterize single shapefile based on its attributes automaticly. Let us suppose that I have shapefile at least 50 attributes. I would like to rasterize this shapefile based on its unigue field attribute in automatic way. It means that every unigue attribute will be stored in single raster band with values of pixels of this attribute. I tried to write an R script, but I was not successfull. Here is part of my code:

library(raster)
library(rgdal)

# setting working directory
vstupni_adresar <- choose.dir(getwd(), "Zadej cestu k pracovnímu adresáři")
prac_adresar <- setwd(vstupni_adresar)

# loading raster data
r <- brick(choose.files(getwd(), "Vyber obrazová data"))

# loading vector data
shp <- shapefile(choose.files(getwd(), "Vyber polygonová data"))

# resolution and extent of raster
cellsize <- 20
ext <- extent(r)
rs <- raster(ext, res=cellsize)

# loop of rasterization
for (i in 3:names(shp)){

# vlastni proces rasterizace    
rasterizace <- rasterize(shp, field=i, rs)
plot(rasterizace)

# name of output raster
outname <- paste("B", i, ".tif", sep="_")

# export raster
rase <- writeRaster(rasterizace, filename=outname, format='GTiff',
                    overwrite=T)
}

But script produces an error. For the matter of speeding up computations I would like to add parallel implementation with foreach package as well.

As for the task itself I does not have to be R, I tried to do this with gdal_rasterize utility which is really fast, but I did not succeed. I prefer R to anything else because it is the thing of which I am the most familiar with.

In another way. For example when I want to rasterize polygon shapefile with attribute called ,,AREA" and save it as a raster file named ,,RASTER" it will look like this:

gdal_rasterize -a AREA -tr 20 20 -of Gtiff shapefile.shp RASTER.tif

But there is a situation when one would like to rasterize the whole attribute table it means to create separate command or batch file with commands above where only parameter -a will differ for the whole attribute table based on name of attribute column. Simplier and much more effective way to do so is to create a script which will loop through list of all present atrributes of polygon shapefile.

  • What went wrong with gdal_rasterize? – user30184 Jun 20 '17 at 17:20
  • Check out fasterize – mdsumner Jun 21 '17 at 6:42

Your Answer

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

Browse other questions tagged or ask your own question.