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

Access Parameter Values from getTerraClim Output for 900 Counties #5

Closed MengyaTao closed 5 years ago

MengyaTao commented 5 years ago

Hi Mike,

I am trying to use function 'getTerraClim', but I am not familiar with the output data structure. Do you know what the commands would be if I want to access parameter values only? Is this something like 'tc$prcp@something'? So I can then be able to save them in .csv in R.

Here is the google drive link to all the counties (.shp file) that I would like to access their parameters from TerraClimate data source: https://drive.google.com/file/d/1y_NN771ch-YS4nSNH6pXes40iELAnR-i/view?usp=sharing

Thank you! Mengya

mikejohnson51 commented 5 years ago

Hi,

From my understanding, you are looking to get an areal average for each of your 891 counties from the TerraClim dataset for all months in 2010. To do this I would loop through the counties using the for each package, here is an example for 'aet'. The trick here is to account for the area of each county yielding a mm/km2 value:

#Define your path 
path = './corn_890_counties.shp'
# Read in your county object
z = read_sf(path)
#Specify a parameter
params = 'aet'

system.time({
    `%dopar%` <- foreach::`%dopar%`
    no_cores  <- parallel::detectCores() - 1
    doParallel::registerDoParallel(no_cores)

    var = foreach::foreach(i = 1:NROW(z),.combine = 'rbind', .packages = 'climateR') %dopar% {
      #assign area in km2                  
      if(param %in% c('tmax', 'tmin')) { area = 1 } else { area = as.numeric(st_area(z[i,]) / 1e6)}

      #get data for county in 2010
      r = getTerraClim(z[i, ],  param = params,  startDate = '2010')

       # calculate mean of each raster layer and append to county FIP code     
      c(as.numeric(z$GEOID[i]),  (cellStats(r[[1]], stat = 'mean', na.rm = TRUE) * area))
     }

    # rename
    colnames(var) = c("GEOID", month.abb)
    row.names(var) = NULL
    # write to csv in the same folder as your shp
    write.csv(var, paste0(dirname(path), "/", params, ".csv"))
})

   user  system elapsed 
 56.238   8.957 114.668 

Hope that helps!

Mike

MengyaTao commented 5 years ago

Thank you so much!!!