2

I would like to plot a raster in WGS84 projection using R based on its attribute value.

Here is the whole procedure

library(maptools)
library(rgdal)
library(raster)

shp <- readOGR("final_result.shp")
cell_size <- 0.012
lon_min <- 10.8; lon_max <- 18.2; lat_min <- 42; lat_max <- 48
ncols <- ((lon_max - lon_min)/cell_size)+1; nrows <- ((lat_max -  lat_min)/cell_size)+1
snap1 <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=0.012, crs="+proj=longlat +datum=WGS84")
snap2 <- raster(nrows=nrows, ncols=ncols, xmn=lon_min, xmx=lon_max, ymn=lat_min, ymx=lat_max, res=0.012, crs="+proj=longlat +datum=WGS84")
snap1v <- shp[shp$snap == 1,]
print(snap1v)
     class       : SpatialPolygonsDataFrame 
     features    : 1 
     extent      : 10.9809, 18.18933, 42.07696, 47.89763  (xmin, xmax, ymin, ymax)
     coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
     variables   : 8
     names       : snap,CO,NH3,NOx,SOx,PM25,PM10,VOC 
     min values  :  1, 0.0510012180001,0,1.25426065201,1.671099876,0.049, 0.0735,0 

snap2v <- shp[shp$snap == 2,]
r_snap1 <- rasterize(snap1v, snap1, fields=c("CO", "NH3", "NOx", "SOx", "PM25", "PM10", "VOC"))

And description of the raster file:

r_snap1
   class       : RasterLayer 
   dimensions  : 500, 617, 308500 (nrow, ncol, ncell)
   resolution  : 0.012, 0.012  (x, y)
   extent      : 10.9809, 18.18933, 42.07696, 47.89763   (xmin, xmax, ymin, ymax)
   coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
   data source : in memory
   names       : layer 
   values      : 1, 1  (min, max)
   attributes  :
       ID snap         CO NH3      NOx    SOx  PM25   PM10 VOC
        1    1 0.05100122   0 1.254261 1.6711 0.049 0.0735   0

How can I plot the attribute value e.g. CO? Just using plot(r_snap1) plots value 1 (ID).

2

How silly of me. Just found the answer:

You only need to add y value

plot(r_snap1, "CO")
0

Here's a ggplot2 approach fwiw, untested code.

## ... prep as per OP
tab <- as.data.frame(r_snap1, xy = TRUE)
library(ggplot2)
ggplot(tab, aes(x, y, fill = CO)) + geom_raster()

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.