4

I've got a shapefile that I would like to convert into a raster, but the trick is that the attribute I want a raster version of is character data. I'm familiar with geoprocessing using GDAL (in the command line) and R, so solutions using either of those two are preferred. However, I have access to QGIS if necessary, and even ArcMap.

I'm using this shapefile. It contains polygons of municipalities in part of Brazil. The attribute NAME_2 is the name of the municipality, and I want a raster version of that.

I know that you can't make a raster with characters, so I created a numeric version of NAME_2 in R by doing

library("rgdal")
munisCrop.shape = readOGR(dsn=getwd(), layer="munisCrop")      # Read shapefile
munisCrop.shape$NAME_2_NUM = as.numeric(munisCrop.shape$NAME_2)  # Make new attribute
writeOGR(obj = munisCrop.shape, dsn = getwd(), layer = "munisCrop", driver = "ESRI Shapefile") # Save new version of shapefile

Values in munisCrop.shape$NAME_2_NUM range from 1 to 787.

I was hoping to then use gdal_rasterize on my command line (GDAL version 1.10.1 on Mac OS X 10.9.2) to make a raster version of the shapefile. Here's the code I was trying (separate lines for clarity):

gdal_rasterize -of GTiff -ot UInt32 -a NAME_2_NUM \
               -tr 0.0002777778 0.0002777778 -ot BYTE \
               -co COMPRESS=LZW -l munisCrop \
               munisCrop.shp munisCrop.tif

This completes easily enough. However, when viewing munisCrop.tif in QGIS, I see that the data range from 0 to 255. This is confirmed in R with

library("raster")
munisCrop.raster = raster("munisCrop.tif")
sort(unique(getValues(munisCrop.raster)))

0 makes sense, as NA values got converted to that. However, the maximum value of 255 is no good; remember, the maximum value of that attribute in the shapefile is 787.

I think I understand why this is happening, but I'd appreciate if somebody could really explain it to me and suggest a workaround.

I've tried various options in gdal_rasterize instead of -ot UInt32 (for example, BYTE), and the result is always the same. I've also tried rasterizing it in R:

munisCrop.blankraster = raster(ext = extent(munisCrop.shp), 
                               resolution = 0.0002777778, 
                               crs = projection(munisCrop.shp) 
                               ) 
munisCrop.Rraster = rasterize(x = munisCrop.shp, 
                              y = munisCrop.blankraster, 
                              field = "NAME_2_NUM" 
                              )

But that seems like it takes forever. gdal_rasterize is way faster (a matter of seconds instead of... indefinite), so I'd love a solution using that.

Thanks in advance!

2

255 is the default NoData value in QGIS.

I am not sure, what exactly the problem is with the way you tried it, but you could use the GDAL Python bindings to do what you want. For instance the following script converts your shp to a polygon based on the attribute NAME_2_NUM.

Import the libraries

import ogr, gdal, osr

Open your shapefile

source_ds = ogr.Open("munisCrop/munisCrop.shp")
source_layer = source_ds.GetLayer()

Create the raster file

pixelWidth = pixelHeight = 0.01 # depending how fine you want your raster
x_min, x_max, y_min, y_max = source_layer.GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)
target_ds = gdal.GetDriverByName('GTiff').Create('yourRaster.tif', cols, rows, 1, gdal.GDT_Byte) 
target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
band = target_ds.GetRasterBand(1)
NoData_value = 999999
band.SetNoDataValue(NoData_value)
band.FlushCache()

Rasterize your sip based on the attribute Name_2_NUM

gdal.RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=NAME_2_NUM"])  

Add a spatial reference

target_dsSRS = osr.SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds.SetProjection(target_dsSRS.ExportToWkt())

enter image description here

0

I have no problem with the speed of the calculation with the function rasterize:

library(raster)
## Set up a raster 
#Extent
ext <-  extent(munisCrop.shape)
#Resolution
xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
r <- raster(ext, ncol=xy[1]/0.1, nrow=xy[2]/0.1)

## Rasterize the shapefile
#you need to define the the value(s) to be transferred (munisCrop.shape$NAME_2_NUM)
# or else the raster values will consist of each unique polygon id
rr <-rasterize(munisCrop.shape, r,field=munisCrop.shape$NAME_2_NUM)
# Plot
plot(rr)
  • Right, as outlined in my post I tried that, but it seemed like it was taking forever. That's why I hoped a GDAL-based solution would be possible and faster. – Sam Rabin Apr 8 '14 at 13:05
  • It didn't took forever. what machine are you using? – Gago-Silva Apr 8 '14 at 16:17

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.