I cropped climate data obtained from Chelsa to obtain a raster for the Netherlands using a country shape file. Plotting both indicates that they are misaligned. The projection is the same.
I followed the introduction on raster on neonscience.org and have read several posts and package documentations in the last two day, as it turns out to be surprisingly difficult to get climate data for GPS data. (Or I am just more incompetent than I have hoped)
setwd(pathtodata)
library(rgdal)
library(raster)
# data for raster http://chelsa-climate.org/downloads/ (adjust file name)
climate <- raster('CHELSAcruts_prec_8_2015_V.1.0.tif')
# data for NUTS shape files https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts#nuts16
bounds <- readOGR(dsn=getwd(),'NUTS_RG_01M_2016_4326_LEVL_0')
bounds <- bounds[bounds$NUTS_ID=='NL',]
climate_nl <- crop(x=climate,y=bounds)
climate_nl@data@values[which(climate_nl@data@values==-32768)] <- NA
plot(climate_nl)
plot(bounds,add=T)
The result would crop the original .tif file to only depict the Netherlands. I would like to then continue with extracting values for GPS locations. Currently, I obtain NAs for locations that should be within the extent.
Like you I got many NAs but I was able to get all but one if I buffered the bounds
extent:
# remotes::install_github("adamhsparks/GSODRdata", build_vignettes = TRUE)
library(GSODRdata)
library(GSODR)
library(dplyr)
library(sf)
library(raster)
climate <- raster("CHELSAcruts_prec_8_2015_V.1.0.tif")
# https://stackoverflow.com/a/49159943/3362993
climate <- reclassify(climate, cbind(-Inf, 0, NA), right=FALSE)
NL <- get_GSOD(years = 2010, country = "Netherlands", max_missing = 5)
bounds <- GSODRdata::CHELSA[,c("STNID", "CHELSA_prec_8_1979-2013")] %>%
left_join(NL, by = "STNID") %>%
dplyr::filter(!is.na(LON)) %>%
distinct(STNID, `CHELSA_prec_8_1979-2013`, .keep_all = TRUE) %>%
dplyr::filter(LON > -33) %>% # fix misclassified county codes
st_as_sf(coords = c("LON", "LAT"), crs = 4326)
climate_nl <- raster::crop(climate, as_Spatial(st_buffer(bounds, 0.1)))
res <- extract(climate_nl, as_Spatial(bounds))
plot(climate_nl)
plot(bounds, add = TRUE, color = "black")
plot(bounds[which(is.na(res)),]$geometry, add = TRUE, color = "red", pch = 19)
legend("topleft", legend = c("NA vals"), pch = 19)