5

I want to plot a shapefile over a raster file in R but I can't make them overlap perfectly: the raster appears to be rotated of few degrees counter-clockwise. Is it a problem with the projection?

enter image description here

Please consider the following example:

library(raster)
library(rgdal) 

# Download from http://biogeo.ucdavis.edu/data/gadm2/shp/ITA_adm.zip
shape_file = "ITA_adm1.shp"
# Download from http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density/data-download
# Setting Geography: Country, Italy; Data Attributes: Grid
pop_density_file ="w001001.adf"

italy_map <- readOGR(dsn = shape_file, layer = "ITA_adm1")
italy_map_dens <- raster(pop_density_file)

colPal <- colorRampPalette(c("white", "red"))( 500 )

par(mar=c(0,0,0,0))
plot(italy_map_dens, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), axes=FALSE, box=FALSE, legend=FALSE, col=colPal)
plot(italy_map, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), border="grey", add=TRUE)
  • Not much help but just a comment: both datasets are in plain wgs84 and work perfectly fine in QGIS. +proj=longlat +datum=WGS84 +no_defs . I don't know much about displaying in R but your screenshot doesn't look like plain WGS84. It looks to me like there's some weird transformation going on. – jareks Nov 21 '13 at 11:14
  • 1
    I would bet that you resized the window after plotting. For some reason the vector gets skewed when you resize the plotting device. Try plotting to a difference graphic device (i.e., jpeg, pdf) and see if the problem persist. – Jeffrey Evans Nov 21 '13 at 15:54
6

I could not reproduce your problem in R (similarly to jareks, in QGIS). My output map was ok, apparently.

This is data from year 2000.

enter image description here

This is your code (adapted) to reproduce all the three years (1990,1995,2000) of population density available on the very interesting link you provided. They seem all fine.

library(raster)
library(rgdal) 

italy_map <- readOGR(dsn = "C:\\...\\ITA_adm", layer = "ITA_adm1")

#ds90g  population densities in 1990, unadjusted, persons per square km
pop_density_file_1990 ="C:\\...\\ita_gpwv3_pdens_wrk_25\\itadens\\itads90ag\\w001001.adf"
#ds95g  population densities in 1995, unadjusted, persons per square km
pop_density_file_1995 ="C:\\...\\ita_gpwv3_pdens_wrk_25\\itadens\\itads95ag\\w001001.adf"
#ds00g  population densities in 2000, unadjusted, persons per square km
pop_density_file_2000 ="C:\\...\\ita_gpwv3_pdens_wrk_25\\itadens\\itads00ag\\w001001.adf"

italy_map_dens_1990 <- raster(pop_density_file_1990)
italy_map_dens_1995 <- raster(pop_density_file_1995)
italy_map_dens_2000 <- raster(pop_density_file_2000)

colPal <- colorRampPalette(c("white", "red"))( 500 )

par(mfrow=c(1,3))

par(mar=c(0,0,0,0))
plot(italy_map_dens_1990, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), axes=FALSE,   box=FALSE, legend=FALSE, col=colPal)
plot(italy_map, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), border="grey", add=TRUE)

par(mar=c(0,0,0,0))
plot(italy_map_dens_1995, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), axes=FALSE,  box=FALSE, legend=FALSE, col=colPal)
plot(italy_map, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), border="grey", add=TRUE)

par(mar=c(0,0,0,0))
plot(italy_map_dens_2000, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), axes=FALSE, box=FALSE, legend=FALSE, col=colPal)
plot(italy_map, xlim = c(6.70, 18.32), ylim = c(35.2, 47.6), border="grey", add=TRUE)

enter image description here

4

Apparently there was a bug in the rgdal package. My problem was solved after I updated it to version 0.8-12.

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.