空间矩形坐标中定义矩形顶点的所有可能组合

Question

由于内存问题，我决定通过矩形从EMODnetWCS下载信息，然后合并它们。

library(EMODnetWCS)

### Script for downloading information from https://github.com/EMODnet/EMODnetWCS
## more specifically wcs information for bathymetry

wcs <- emdn_init_wcs_client(service = "bathymetry")

#Working example
emdn_get_wcs_info(wcs)
emdn_get_coverage_info(wcs, 
                       coverage_ids = "emodnet__mean")

然后你可以像这样下载：

cov <- emdn_get_coverage(wcs, 
                         # coverage_id = "emodnet__mean",
                         # bbox = c(xmin = 0, 
                         #          ymin = 40, 
                         #          xmax = 2, 
                         #          ymax = 42),
                         # nil_values_as_na = TRUE)

我需要从 (lat,lon)= (-15,35) 到 (13,63) 的信息。但信息量太大了，所以我需要分部分来做。

我笨拙的做法是这样循环，将其分成1度^2的矩形：

x=seq(-15,13,1)

length(x)

y=seq(35,63,1)

length(y)

xy <- data.frame(x,y)

my_list <- list()

for (i in 1:(dim(xy)[1]-1)){
  print(i)
  
  xmin=xy[i,1]
  ymin=xy[i,2]
  xmax=xy[i+1,1]
  ymax=xy[i+1,2]
  
  bbox=c(xmin=xmin,ymin=ymin,xmax=xmax,ymax=ymax)
  print(bbox)
  
  cov <- emdn_get_coverage(wcs, 
                           coverage_id = "emodnet__mean",
                           bbox = bbox,
                           nil_values_as_na = TRUE)
  
  rast <- raster(cov)
  my_list[[i]] <- rast
  }

但这只会下载对角线（简化的）

那么，我该怎么做才能为每个矩形定义正确的坐标？

Answer 1

要构建具有所有组合的框架，您可以使用

expand.grid()

:

expand.grid(
  x = seq(-15, 13, 1),
  y = seq(35, 63, 1)
  ) 
#>       x  y
#> 1   -15 35
#> 2   -14 35
#> 3   -13 35
#> 4   -12 35
#> 5   -11 35
# ...

---

虽然您已经将

sf

作为

EMODnetWCS

要求之一，但您也可以使用

sf::st_make_grid()

来构建网格单元并使用

sf::st_bbox()

来获取 bboxes :

library(EMODnetWCS)
#> Loading ISO 19139 XML schemas...
#> Loading ISO 19115 codelists...
library(purrr)
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is FALSE
library(terra)
#> terra 1.7.83

# make grid
grid_sfc <- 
  st_bbox(c(xmin = -15, ymin = 35, xmax = 13, ymax = 63)) |> 
  st_as_sfc() |> 
  st_make_grid(cellsize = 1)
grid_sfc
#> Geometry set for 784 features 
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -15 ymin: 35 xmax: 13 ymax: 63
#> CRS:           NA
#> First 5 geometries:
#> POLYGON ((-15 35, -14 35, -14 36, -15 36, -15 35))
#> POLYGON ((-14 35, -13 35, -13 36, -14 36, -14 35))
#> POLYGON ((-13 35, -12 35, -12 36, -13 36, -13 35))
#> POLYGON ((-12 35, -11 35, -11 36, -12 36, -12 35))
#> POLYGON ((-11 35, -10 35, -10 36, -11 36, -11 35))

# check bboxes for first 3 cells;
# plot grid, first 3 cells highlighted
grid_sfc[1:3] |> map(st_bbox)
#> [[1]]
#> xmin ymin xmax ymax 
#>  -15   35  -14   36 
#> 
#> [[2]]
#> xmin ymin xmax ymax 
#>  -14   35  -13   36 
#> 
#> [[3]]
#> xmin ymin xmax ymax 
#>  -13   35  -12   36
plot(grid_sfc, col = c("white", "red")[1+(seq(grid_sfc) <= 3)], axes = TRUE)

wcs <- emdn_init_wcs_client(service = "bathymetry")
#> ✔ WCS client created succesfully
#> ℹ Service: <https://ows.emodnet-bathymetry.eu/wcs>
#> ℹ Service: "2.0.1"
cov_lst <-
  # download first 3 cells 
  grid_sfc[1:3] |> 
  imap(
    \(grid_cell, idx) emdn_get_coverage(
      wcs = wcs, 
      coverage_id = "emodnet__mean", 
      bbox = st_bbox(grid_cell), 
      filename = sprintf("%.3d.tif", idx), 
      nil_values_as_na = TRUE
    )
  )
#> ── Downloading coverage "emodnet__mean" ────────────────────────────────────────
#> No encoding supplied: defaulting to UTF-8.
#> <GMLEnvelope>
#> ....|-- lowerCorner: 35 -15
#> ....|-- upperCorner: 36 -14Start tag expected, '<' not found
#> 
#> ✔ Coverage "emodnet__mean" downloaded succesfully as a
#>        terra <SpatRaster>
#> ✔ nil values NA converted to NA on all bands.
#> ── Downloading coverage "emodnet__mean" ────────────────────────────────────────
#> <GMLEnvelope>
#> ....|-- lowerCorner: 35 -14
#> ....|-- upperCorner: 36 -13Start tag expected, '<' not found
#> 
#> ✔ Coverage "emodnet__mean" downloaded succesfully as a
#>        terra <SpatRaster>
#> ✔ nil values NA converted to NA on all bands.
#> ── Downloading coverage "emodnet__mean" ────────────────────────────────────────
#> <GMLEnvelope>
#> ....|-- lowerCorner: 35 -13
#> ....|-- upperCorner: 36 -12Start tag expected, '<' not found
#> 
#> ✔ Coverage "emodnet__mean" downloaded succesfully as a
#>        terra <SpatRaster>
#> ✔ nil values NA converted to NA on all bands.
# list of resulting SpatRasters
cov_lst  
#> [[1]]
#> class       : SpatRaster 
#> dimensions  : 960, 960, 1  (nrow, ncol, nlyr)
#> resolution  : 0.001041667, 0.001041667  (x, y)
#> extent      : -15, -14, 35, 36  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> varname     : 001 
#> name        :       001 
#> min value   : -4754.600 
#> max value   : -1791.212 
#> 
#> [[2]]
#> class       : SpatRaster 
#> dimensions  : 960, 960, 1  (nrow, ncol, nlyr)
#> resolution  : 0.001041667, 0.001041667  (x, y)
#> extent      : -14, -13, 35, 36  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> varname     : 002 
#> name        :        002 
#> min value   : -5064.2476 
#> max value   :  -681.6658 
#> 
#> [[3]]
#> class       : SpatRaster 
#> dimensions  : 960, 960, 1  (nrow, ncol, nlyr)
#> resolution  : 0.001041667, 0.001041667  (x, y)
#> extent      : -13, -12, 35, 36  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source(s)   : memory
#> varname     : 003 
#> name        :         003 
#> min value   : -4951.20947 
#> max value   :   -98.78173

# stored files
fs::dir_info(glob = "0*.tif")[,1:5]
#> # A tibble: 3 × 5
#>   path       type         size permissions modification_time  
#>   <fs::path> <fct> <fs::bytes> <fs::perms> <dttm>             
#> 1 001.tif    file           4M rw-         2024-10-16 16:29:45
#> 2 002.tif    file           4M rw-         2024-10-16 16:29:49
#> 3 003.tif    file           4M rw-         2024-10-16 16:29:54

^{创建于 2024 年 10 月 16 日，使用 reprex v2.1.1}