1

有什么方法可以用来解析一个国家的 shapefile 并使用 R 在该国家/地区下载 MODIS 产品数据？

我尝试了使用MODIStsp包（https://docs.ropensci.org/MODIStsp/）以及MODISTools包（ https://docs.ropensci.org/MODISTools/articles/modistools-vignette.html ）的不同方法，它们都只允许我下载特定站点的 MODIS 产品数据，但不能下载国家/地区。

4

2 回答 2

2

这是一个如何实现这一目标的示例。

首先，下载你需要的 MODIS 数据，在这个例子中我使用的是 MCD12Q1.006

begin_yearend_year格式为：Year.Month.Days 。shape_file是您正在使用的 shapefile，大概 shapefile 的范围是您所追求的国家/地区。不过，我只是根据您提供的最少信息而离开。

library(MODIS)

tifs <- runGdal(product = "MCD12Q1", collection = "006", SDSstring = "01", 
                extent = shape_file %>% st_buffer(dist = 10000), 
                begin = begin_year, end = end_year, 
                outDirPath = "data", job = "modis",
                MODISserverOrder = "LPDAAC") %>% 
  pluck("MCD12Q1.006") %>% 
  unlist()

# rename tifs to have more descriptive names
new_names <- format(as.Date(names(tifs)), "%Y") %>% 
  sprintf("modis_mcd12q1_umd_%s.tif", .) %>% 
  file.path(dirname(tifs), .)
file.rename(tifs, new_names)

landcover <- list.files("data/modis", "^modis_mcd12q1_umd", 
                        full.names = TRUE) %>% 
  stack()
# label layers with year
landcover <- names(landcover) %>% 
  str_extract("(?<=modis_mcd12q1_umd_)[0-9]{4}") %>% 
  paste0("y", .) %>% 
  setNames(landcover, .)

此外，如果您需要特定的像元大小，则可以按照此过程获得 5x5 modis 像元大小。

neighborhood_radius <- 5 * ceiling(max(res(landcover))) / 2

agg_factor <- round(2 * neighborhood_radius / res(landcover))
r <- raster(landcover) %>% 
  aggregate(agg_factor) 
r <- shape_file %>% 
  st_transform(crs = projection(r)) %>% 
  rasterize(r, field = 1) %>% 
  # remove any empty cells at edges
  trim()
于 2021-03-03T12:07:22.630 回答
1

这是一个使用 MODISTools 自动下载该国家/地区的正确图块的示例。

首先让我们生成一个国家的多边形来演示（以卢森堡为例）：

library(maptools)
library(sf)
data(wrld_simpl)
world = st_as_sf(wrld_simpl)
lux = world[world$NAME=='Luxembourg',]

现在我们找到国家的位置（质心）和大小：

#find centroid of polygon in long-lat decimal degrees
lux.cent = st_centroid(lux)

#find width and height of country in km
lux.proj = st_transform(lux, 
           "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=km +no_defs")
lux.km_lr = diff(st_bbox(lux.proj)[c(1,3)])
lux.km_ab = diff(st_bbox(lux.proj)[c(2,4)])

使用此信息，我们可以下载正确的 Modis 数据（以叶面积索引 lai 为例）：

#download the MODIS tiles for the area we defined
library(MODISTools)
lux_lai <- mt_subset(product = "MOD15A2H",
                          lat = lux.cent$LAT, lon =  lux.cent$LON,
                          band = "Lai_500m",
                          start = "2004-01-01", end = "2004-01-01",
                          km_lr = lux.km_lr, km_ab = lux.km_ab,
                          site_name = "Luxembourg",
                          internal = TRUE, progress = TRUE)

# convert to a spatial raster
lux.rast = mt_to_raster(df = lux_lai, reproject = TRUE)
lux.rast = raster::mask(lux.rast, lux)
plot(lux.rast)
plot(st_geometry(lux),add=T)

在此处输入图像描述

于 2021-03-06T16:18:51.697 回答