2

I am having trouble converting a polygon to a raster. The desired raster is created, however all its values are equal to NA.

require(rgdal)
require(raster)

#download a shapefile with 2010 census data
tmp_dl <- tempfile()
download.file("http://files.hawaii.gov/dbedt/op/gis/data/blkgrp10.shp.zip", tmp_dl)
unzip(tmp_dl, exdir=tempdir())
HIshp <- readOGR(tempdir(), "blkgrp10")

#add a new attribute/field with population density of each polygon
p_areas <- sapply(HIshp@polygons, function(i)i@area)
p_pops <- HIshp@data$POP10
HIshp@data[,3] = p_pops/p_areas
names(HIshp@data)[3] <- "POPDENS"

#here is the raster format that I require for the result
new_ras <- raster(nrow = 1520, ncol = 2288)
extent(new_ras) <- c(-159.816, -154.668, 18.849, 22.269)
crs(new_ras) <- "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"

rasterize(HIshp, new_ras, field="POPDENS")

And here is the sad result:

> rasterize(HIshp, new_ras)
class       : RasterLayer 
dimensions  : 1520, 2288, 3477760  (nrow, ncol, ncell)
resolution  : 0.00225, 0.00225  (x, y)
extent      : -159.816, -154.668, 18.849, 22.269  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=4 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : NA, NA  (min, max)

How can I make the values convert properly?

  • 1
    It is projected shapefile NAD_1983_UTM_Zone_4N. Now overlaps of course with required raster extent. – FelixIP Feb 9 '15 at 1:27
2

First, reproject the vector data to WGS84 (Lat/Long degrees):

HI_WGS84 <- spTransform(HIshp, CRS("+proj=longlat +elips=WGS84"))

Then rasterize a new result:

new_ras <- raster(nrow=1520, ncol=2288)
crs(new_ras) <- crs(HI_WGS84)
extent(new_ras) <- c(-159.816, -154.668, 18.849, 22.269)
HI_popdens <- rasterize(HI_WGS84, new_ras, field="POPDENS")
plot(HI_popdens)

plot(HI_popdens)

  • Unfortunately I require the output to be in a different projection from the shapefile. Plus one for the helpful code and images. – J. Win. Feb 12 '15 at 0:34
  • Then use spTransform; see simplified answer. – Mike T Feb 12 '15 at 3:18
2

Well, you have two potential issues:

1) "HIshp" is in a "UTM Zone 4" projection and the raster you defined is in a geographic projection (lat/long). You can use "spTransform" to reproject "Hishp" so that it aligns with "new_ras".

2) You need to define the attribute in "HIshp" that you want to represent the raster values (ie., rasterize(HIshp, new_ras, field="POPDENS").

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.