1

I posted a question (please find the link below) regarding raster plots at stackoverflow but could not get the solution. Does anyone here has the solution? https://stackoverflow.com/questions/63209707/how-to-perform-point-to-raster-operation-and-plot-it-with-raster-of-different-sp

The question :

I have following dataframe with three columns namely, "Latitude", "Longitude", and "value", along with a raster file (as attached with name "ras.tif') and area shape file (as attached with name "okl_shape.shp"):

Dataframe

df <- data.frame(Latitude=c(34.80833, 34.79851, 34.91418, 36.07204, 36.80253, 35.96305, 
36.69256, 36.83129, 35.5915, 34.8497, 34.60896, 35.65282, 36.74813, 35.54615, 33.92075, 
35.54848, 35.20494, 36.26353, 36.84053, 36.60183), Longitude=c(-98.02325, -96.66909, -98.29216, 
-99.90308, -100.53012, -95.86621, -102.49713, -99.64101, -99.27059, -97.0033, -96.33309, 
-96.80407, -98.36274, -99.7279, -96.32027, -98.03654, -99.80344, -98.49766, -96.42777, -101.6013 
), value=c(0.84629845, NA, 0.916147287, 0.735707364, 0.432443798, 0.959728682, 0.123924419, 
0.849589147, 0.307998062, 0.932215116, 0.939897287, 0.915825581, NA, 0.273945736, NA, 
0.705023256, 0.447494186, 0.840686047, 0.901098837, 0.202523256))

Attached files

Here is the shape file of my area: file Here is the raster file: file

Problem Description

I have to make a plot comparing point data (i.e. data frame described above) with attached raster (i.e. "ras.tif") and the shape file should overlay these plots. See the picture below for reference (In this figure, the left panel plots the points values as defined in the dataframe, right panel plots the raster, and the polygon represent the area shape file): enter image description here

My code

I tried solving this with "spplot" but couldn't get the desired plot (something similar to above plot). Also, I don't know how can I insert the shape file in my plot. Below is the code I tried so far:

library(raster)
library(sf)
library(maptools)
meso.shape <- st_read('okl_shape.shp')
ras1 <- raster('ras.tif')
df <- data.frame(Latitude=c(34.80833, 34.79851, 34.91418, 36.07204, 36.80253, 35.96305, 
36.69256, 36.83129, 35.5915, 34.8497, 34.60896, 35.65282, 36.74813, 35.54615, 33.92075, 
35.54848, 35.20494, 36.26353, 36.84053, 36.60183
), Longitude=c(-98.02325, -96.66909, -98.29216, -99.90308, -100.53012, -95.86621, -102.49713, 
-99.64101, -99.27059, -97.0033, -96.33309, -96.80407, -98.36274, -99.7279, -96.32027, -98.03654, 
-99.80344, -98.49766, -96.42777, -101.6013), value=c(0.84629845, NA, 0.916147287, 0.735707364, 
0.432443798, 0.959728682, 0.123924419, 0.849589147, 0.307998062, 0.932215116, 0.939897287, 
0.915825581, NA, 0.273945736, NA, 0.705023256, 0.447494186, 0.840686047, 0.901098837, 
0.202523256))
coordinates(df) <- ~Longitude+Latitude
x <- vect2rast(df, fname = names(df)[1])
x <- raster(x)
extent(x) <- extent(ras1)
crs(x) <- crs(ras1)
x <- projectRaster(x,ras1)
s <- stack(ras1,x)
sp <- as(s, 'SpatialPolygons')
p <- spplot(s, names.attr = c("Simulated FWI","Observed FWI"))
show(p)
1
  • There is no method for supporting raster class objects in spplot. You can coerce to a sp object using as (x, "SpatialPixelsDataFrame") which can then be plotted in the lattice plot engine. If the data is in the same projection just plot your raster plot(ras1 ) and then overlay the points using points(meso.shp, pch=19, cex=0.75) Just like plot', the points` function will take color arguments for symbology. You can do the same thing by calling plot again then including the add=TRUE argument. Sep 10 '20 at 22:37
1

You are very close, although mixing code from spand sflibraries; you only attached the shp file, so I couldn't download your file, but made a dummy linestring object from your points so you can see how to plot a line on a raster. It is easier to use lower functions (par(mfrow = c(2,1)) to plot the line on both rasters of the stack

library(raster)
library(sf)
library(dplyr) # to make sf flow

ras1 <- raster('ras.tif')

df <- data.frame(Latitude=c(34.80833, 34.79851, 34.91418, 36.07204, 36.80253, 35.96305, 
                            36.69256, 36.83129, 35.5915, 34.8497, 34.60896, 35.65282, 36.74813, 35.54615, 33.92075, 
                            35.54848, 35.20494, 36.26353, 36.84053, 36.60183), Longitude=c(-98.02325, -96.66909, -98.29216, 
                                                                                           -99.90308, -100.53012, -95.86621, -102.49713, -99.64101, -99.27059, -97.0033, -96.33309, 
                                                                                           -96.80407, -98.36274, -99.7279, -96.32027, -98.03654, -99.80344, -98.49766, -96.42777, -101.6013 
                            ), value=c(0.84629845, NA, 0.916147287, 0.735707364, 0.432443798, 0.959728682, 0.123924419, 
                                       0.849589147, 0.307998062, 0.932215116, 0.939897287, 0.915825581, NA, 0.273945736, NA, 
                                       0.705023256, 0.447494186, 0.840686047, 0.901098837, 0.202523256))

df = st_as_sf(df, coords = c("Longitude", "Latitude")) %>% st_set_crs(4326)
df = df %>% filter(!is.na(value))
a_line = df %>% summarise() %>% st_cast("LINESTRING") # dummy linestring
x = raster::rasterize(df, 
                      raster::raster(df, resolution = 0.25), 
                      df$value)

ras1 = crop(ras1, extent(x), snap = "near")

ras1 = resample(ras1, x)
s <- stack(ras1,x)
plot(s)
par(mfrow = c(1,2))
plot(x, col = RColorBrewer::brewer.pal(8, "Spectral"))
plot(a_line, add = T)
plot(ras1, col = RColorBrewer::brewer.pal(8, "Spectral"))
plot(a_line, add = T)

enter image description here

1
  • What is the role of "resolution" in rasterize function? I think rasterize transfers points (or other spatial data) values to raster cells which may lead to shifting of some of the points based on point to point distance. I want to plot the points to be centered exactly at their respective lon-lat. My plot should look like the plot I attached in the question (like those circular points in the attached image).
    – raghav
    Aug 12 '20 at 21:01

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.