How to crop gridded data with shapefile

[Ellenz]

The package is installed, but the same problem appeared again and I could not make the 'mask' object
dist is assumed to be in decimal degrees (arc_degrees).
although coordinates are longitude/latitude, st_difference assumes that they are planar
Error in UseMethod("st_as_sf") :
no applicable method for 'st_as_sf' applied to an object of class "c('sfc_POLYGON', 'sfc')"
In addition: Warning message:
In st_buffer.sfc(., 5) :
st_buffer does not correctly buffer longitude/latitude data

[jlacko]

The warnings are okay, but the failure of st_as_sf signals a problem.

What versions of geos, gdal and proj libraries are reported when you load {sf} package?

[Ellenz]

This may give you the information:
Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
Warning message:
package ‘sf’ was built under R version 3.4.4

[jlacko]

Weird... Try this, it should work as well - it has the offending st_as_sf() eliminated (the mask does not need to be of sf class, sfc is good enough for masking).

library(tidyverse)
library(sf)

borders <- readRDS(url("https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/gadm36_IRN_0_sf.rds"))

DFgrid <- data.frame(
  lon=c(59.25, 59.75, 60.25, 60.75, 61.25, 61.75, 57.25, 57.75, 58.25, 58.75, 59.25, 59.75, 60.25, 60.75, 61.25, 61.75, 62.25, 57.25, 57.75, 58.25, 58.75, 59.25, 59.75, 60.25, 60.75, 61.25, 61.75, 62.25, 62.75, 53.25),
  lat=c(25.25, 25.25, 25.25, 25.25, 25.25, 25.25, 25.75, 25.75, 25.75, 25.75, 25.75, 25.75, 25.75, 25.75, 25.75, 25.75, 25.75, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.25, 26.75),
  value= rnorm(30)
)


mask <- borders %>% 
  st_bbox() %>% 
  st_as_sfc() %>% 
  st_buffer(5) %>% 
  st_difference(borders)


my_year <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())
my_fill <- scale_fill_distiller('',palette='Spectral',
                                breaks = c(-2,0,2))

fig = ggplot() + geom_tile(data = DFgrid, aes(y=lat, x=lon, fill=value)) + my_year + my_fill +
  geom_sf(data = mask, fill = "white", color = NA) +
  coord_sf(xlim=c(44.05, 63.35), ylim=c(25.05, 39.75)) +
  borders('world', xlim=c(44.05, 63.35), ylim=c(25.05, 39.75), colour='black')+
  xlab('Longitude (°E)')+ ylab('Latitude (°N)')+
  theme(panel.grid=element_blank(),
        panel.border = element_rect(colour = "black", fill=NA,size=.7))

print(fig)

[Ellenz]

Thanks a lot, it worked perfectly.
One more question is, the cropped data may be different from (has a smaller range than) the whole data, so is it possible to adapt the color bar to the cropped data. And the crop here uses the second option you gave above.

[jlacko]

This might be a problem; if you avoid cropping on the data level and "just" mask the output at presentation level all the data features (such as range of values) are retained.

Your best option I can think of is:

• first crop at data level via raster::mask() first, to see the range of values
• then do the masking at ggplot2 level, setting the range of legend manually to the values found in step one