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

Many Sites from TerraClim Normals #29

Closed mikejohnson51 closed 3 years ago

mikejohnson51 commented 3 years ago

From an email: I'm trying to use some combination of 3 packages (ncdf4, climateR, raster -- or others if you recommend them!) to extract data for ~40k US sites, for a number of terraclimate normal variables. 


I will soon add the TerraClim normals to climateR. Until then, I would suggest using the raw URLs as follows. Depending of your proj/gdal you may get warnings about CRS elements, this are ok.

library(raster)
library(AOI)
library(sf)
library(dplyr)

random_pts = AOI::aoi_get(state = "conus") %>% 
            st_sample(40000) %>% 
            st_cast("POINT") %>% 
            st_as_sf() %>% 
            mutate(ID = 1:n())

# Assuming you already have your points, start here:

# Bounding Box 
bb = st_as_sf(st_as_sfc(st_bbox(random_pts)))

rb <- brick("http://thredds.northwestknowledge.net:8080/thredds/dodsC/TERRACLIMATE_ALL/summaries/TerraClimate19611990_tmin.nc")

us_terra = crop(rb, bb)

ext = data.frame(ID = random_pts$ID, extract(us_terra, st_as_sf(random_pts)))

ext[1:2, 1:4]
#>   ID X1961.01.01 X1961.02.01 X1961.03.01
#> 1  1       -12.9       -10.1        -6.2
#> 2  2       -20.1       -16.7        -9.5

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

mikejohnson51 commented 3 years ago

Here is a follow up using the recently added getTerraClimNormals function. This is a beta function for now so if you see anything suspious please file an issue!

Thanks

library(raster)
library(AOI)
library(sf)
library(climateR)
library(dplyr)

# 40K random points in CONUS
random_pts = AOI::aoi_get(state = "conus") %>% 
  st_sample(40000) %>% 
  st_cast("POINT") %>% 
  st_as_sf() %>% 
  mutate(ID = 1:n())

# Assuming you already have your points, start here:
# Pass your points, ideal parameters, period of record and desired months...

system.time({
  ext  = getTerraClimNormals(random_pts, param = "tmin", period = '19611990', month = 1:12)
  ext2 = extract_sites(r = ext, pts = random_pts, id = "ID")
})
#>    user  system elapsed 
#>   6.493   1.435  18.319

plot(ext2$`19611990_tmin`$site_2)

In total we downloaded 816,960 data-points to extract the desired 480,000 data-points in less then 20 seconds 👍

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