mikejohnson51 / climateR

An R 📦 for getting point and gridded climate data by AOI
https://mikejohnson51.github.io/climateR/
MIT License
168 stars 40 forks source link

getGridMET returning inconsistent data #31

Closed sgwinder closed 3 years ago

sgwinder commented 3 years ago

Running getGridMET seems to be returning faulty data, particularly towards the end of long time series. The data returned is also inconsistent between identical calls.

For example, the following call

start_date <- "2010-01-01"
stop_date <- "2021-03-01"
site_bbox <- st_bbox(sites)

weather_rast <- getGridMET(AOI = site_bbox, 
                             param = c("prcp", "tmin", "tmax"), 
                             startDate = start_date, 
                             endDate = stop_date)

Runs without any errors (though it does return the warning In CPL_crs_from_input(x) : GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.).

But, when I try examining the data from a particular date using rasterVis::levelplot(weather_rast$tmax$X2019.01.01)

I get a different result every time I rerun the getGridMET call. Some examples of the plots (all generated using the exact code above) are here:

image image image

Earlier dates in the time series do show believable weather data.

Any insights are welcome!

mikejohnson51 commented 3 years ago

Hi @sgwinder thanks for the note, in case you haven't please update your climateR version.

Otherwise I am happy to work through this with you but currently can't reproduce the errors you are seeing, please see below:

library(AOI)
library(climateR)
library(sf)
library(rasterVis)
library(dplyr)

# Random region and points
region <- AOI::aoi_get(state = "CA", county = "Santa Barbara")
sites = sf::st_sample(region, 100)
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

# pass points directly to getGridMet
weather_rast <- climateR::getGridMET(AOI = sites, 
                           param = c("prcp", "tmin", "tmax"), 
                           startDate = "2010-01-01", 
                           endDate = NULL)

# look at a plot
rasterVis::levelplot(weather_rast$tmin)

dplyr::filter(param_meta$gridmet, common.name == "tmax")
#>   common.name call               description units
#> 1        tmax tmmx daily_maximum_temperature  degK

# units are in degree kelvin so seem resonalbe...

# Lets make sure the same results are returned? Test the same call 10 times
test = list()
for(i in 1:10){
  test[[i]] <- climateR::getGridMET(AOI = sites, 
                             param = 'tmax', 
                             startDate = "2010-01-01", 
                             endDate = NULL)[[1]]
}

g = expand.grid(1:10, 1:10)
head(g)
#>   Var1 Var2
#> 1    1    1
#> 2    2    1
#> 3    3    1
#> 4    4    1
#> 5    5    1
#> 6    6    1

check = list()
for(i in 1:nrow(g)){
  check[[i]] = identical(test[[g[i,1]]], test[[g[i,2]]])
}

# Are the TRUES equal to the comparision?
sum(unlist(check)) == nrow(g)
#> [1] TRUE

Created on 2021-03-12 by the reprex package (v0.3.0)

sgwinder commented 3 years ago

Hi @mikejohnson51 , thanks for the quick response! I just updated, but am still getting the same issue.

Can you try pulling a long time series? My early dates are fine, but at some point the data switch to being not interpretable. In this latest attempt, you can see this happens mid-way through getting tmax data for 2016-05-30 (my timeseries begins at 2010-01-01 and attempts to go until 2021-03-01). image

Other variables fail at different times - tmin switches to every value being 210 on 2012-01-20, while prcp actually seems to have pulled good data all the way through the end of the time series on 2021-03-01 (this time, previous attempts haven't gotten me complete time series of prcp either).

Thanks!

mikejohnson51 commented 3 years ago

Hi @sgwinder, would you mind sharing what you are using as sites (save with dput?) ?

I am still getting expected results from what you describe:

library(AOI)
library(climateR)
library(sf)

# Random region and points
region <- AOI::aoi_get(state = "CA", county = "Santa Barbara")
sites = sf::st_sample(region, 1)

start_date <- "2010-01-01"
stop_date <- "2021-03-01"

weather_rast <- getGridMET(AOI = sites, 
                           param = c("prcp", "tmin", "tmax"), 
                           startDate = start_date, 
                           endDate = stop_date)

plot(weather_rast$date, weather_rast$tmin)

plot(weather_rast$date, weather_rast$tmax)

plot(weather_rast$date, weather_rast$prcp)

Created on 2021-03-12 by the reprex package (v0.3.0)

sgwinder commented 3 years ago

Strange. I am using a bounding box as my geographic input - I'm attaching the result of dput here. But in case that doesn't work, the actual site file is also attached (compressed). From the raw sites, I'm doing:

library(sf)
library(tidyverse)
library(climateR)

sites <- read_sf("allsites-wwa_boundaries_20200506.geojson")
site_bbox <- st_bbox(sites)

Then using site_bbox as my AOI. site_bbox_dput.txt

allsites-wwa_boundaries_20200506.zip

mikejohnson51 commented 3 years ago

I have been messing with this some today and am getting consistent reasonable results for the area. levelplot from rasterVis on occasion hung up on rendering the plot, but that doesn't explain what you're seeing.

Maybe try to pass your regions directly?

library(raster)
library(climateR)
library(sf)

region = read_sf('/Users/mikejohnson/Downloads/allsites-wwa_boundaries_20200506.geojson')

weather_rast <- getGridMET(AOI = region, 
                           param = c("tmax", "tmin", "prcp"), 
                           startDate = "2010-01-01", 
                           endDate   = "2021-03-01")

plot(weather_rast$tmax$X2019.01.01)

plot(weather_rast$tmin$X2019.01.01)

plot(weather_rast$prcp$X2019.01.01)

Created on 2021-03-13 by the reprex package (v0.3.0)

sgwinder commented 3 years ago

Hi @mikejohnson51 , thanks for your suggestions - passing the regions directly didn't fix it for me. I've asked someone on my team to see if she can reproduce my issue and I'll follow up if she is able to.

Thanks for the quick feedback!

mikejohnson51 commented 3 years ago

@sgwinder of course, I just closed this to avoid worrying others since I can't reproduce it. If the error persists please reopen or start a new issue!