3

我有一个纬度数据框，我想从“NASA/NEX-GDDP”数据集中提取每日降雨数据，该数据集是每日时间分辨率和 0.25 度空间分辨率。我的思考过程是

  1. 在我的本地驱动器上下载每日栅格。所以对于给定的年份，将有 365 个栅格
  2. 为我的经纬度提取每个每日栅格的值
  3. 删除每日栅格以节省磁盘空间

我需要这样做多年。上述步骤非常耗时，尤其是步骤 1。有什么方法可以直接提取每日数据，而不必下载栅格。以下是我目前的工作方式

  library(rgee)
  ee_Initialize(email = 'xxxx@gmail.com')
  
  # define a bounding box within which my lat lon lies 
  polygon <- ee$Geometry$Polygon(
    list(
      c(-76.12, -36.69),
      c(-15.42,-36.69),
      c(-15.42, 6.71),
      c(-76.12, 6.71),
      c(-76.12, -36.69)
    ))
  
  latlon <- structure(list(lat = c(-11.125, -11.125, -10.875, -10.875, -10.875, -10.875, -10.875, -10.875), 
                           lon = c(-37.375, -37.125, -70.625, -70.375, -70.125, -69.875, -69.625, -69.375)), 
                      class = "data.frame", row.names = c(NA, -8L))
  
  model_ref <- 'ACCESS1-0' 
  rcp_ref <- 'historical'
  year_ref <- 2005
  
  Imagecollection_pr <-
    ee$ImageCollection("NASA/NEX-GDDP")$ # read the data
    select('pr')$ # select rainfall
    filter(ee$Filter$calendarRange(year_ref, year_ref, "year"))$ # filter the year
    filter(ee$Filter$eq("model", model_ref))$ # filter the model
    filter(ee$Filter$eq("scenario", rcp_ref)) # filter the scenario

  # Downloading as raster from Earth Engine to Local. This downloads 365 images for 2005
  ee_imagecollection_to_local(
    ic = Imagecollection_pr,
    region = polygon,
    dsn = paste0('pr'),
    via = 'getInfo')
  
  # extract values for the given latlon
  library(raster)
  library(rgdal)
  
  ll <- latlon
  
  coordinates(ll) <- ~ lon + lat
  proj4string(ll) <- "+proj=longlat +datum=WGS84 +no_defs"
  
  raster_ls <- list.files(file.path(getwd()), pattern = '.tif')
  
  temp_ls <- list()
  
  for(rr in seq_along(raster_ls)){

    raster_temp <- raster_ls[rr]
    raster_doy <- raster(raster_temp)
    
    crs(raster_doy) <- "+proj=longlat +datum=WGS84 +no_defs"
    
    ymd <- 
      as.numeric(gsub(model_ref,"",raster_temp) %>% 
                   gsub('pr',"",.) %>% 
                   gsub(rcp_ref,"",.) %>%
                   gsub(".tif","",.) %>%
                   gsub("___","",.)) 
    
    latlon$val <- extract(raster_doy, ll)
    temp_ls[[rr]] <- latlon %>% dplyr::mutate(ymd = ymd)
  }  
  
  file.remove(raster_ls)
  results <- do.call('rbind', temp_ls)
  head(results)      
  
  lat     lon          val      ymd
  1 -11.125 -37.375 1.677004e-05 20050101
  2 -11.125 -37.125 1.582177e-05 20050101
  3 -10.875 -70.625 4.219083e-05 20050101
  4 -10.875 -70.375 3.825104e-05 20050101
  5 -10.875 -70.125 3.438132e-05 20050101
  6 -10.875 -69.875 2.946513e-05 20050101
4

0 回答 0