r - 如何在全局文件中找到像素(单元格)的 4 个坐标(lat-long)？

Question

我正在使用可以从这里下载的气候变量:

  ftp://sidads.colorado.edu/pub/DATASETS/nsidc0301_amsre_ease_grid_tbs/global/ 

此文件是一个二进制(矩阵)文件，具有 586 行和 1383 列(全局 map )。 我想知道一个像素(单元格)的 4 个坐标(经纬度)`。

关于文件的更多信息:

These data are provided in EASE-Grid projections global cylindricalat 25 km resolution, are two-
      byte 
    Spatial Coordinates:
    N: 90°     S: -90°     E: 180°     W: -180°    

使用栅格包并将数据转换为栅格对象:

 file<- readBin("ID2r1-AMSRE-ML2010001A.v03.06H", integer(), size=2,  n=586*1383, signed=T)
 m = matrix(data=file,ncol=1383,nrow=586,byrow=TRUE)
 r = raster(m, xmn=-180, xmx=180, ymn=-90, ymx=90)

现在该文件是一个正确的空间引用对象，但如果没有使用的圆柱投影的完整规范，您将无法返回经纬度坐标。

这里有更多信息 http://nsidc.org/data/ease/tools.html包括指向某些网格的链接，这些网格具有该网格系统的网格单元的经纬度:

ftp://sidads.colorado.edu/pub/tools/easegrid/lowres_latlon/

      MLLATLSB  "latitude"
      MLLonLSB  "longitude"

例如，我们可以为数据网格中的单元格创建一个纬度栅格:

> lat <- readBin("MLLATLSB",integer(), size=4,  n=586*1383, endian="little")/100000
> latm = matrix(data=lat,ncol=1383,nrow=586,byrow=TRUE)
> latr = raster(latm, xmn=-180, xmx=180, ymn=-90, ymx=90)

然后latr[450,123]是我数据中单元格[450,123]的纬度。对经度重复 MLLONLSB。

但是这还不够(每个像素一个纬度和一个经度)，因为我想与基于地面的测量结果进行比较，所以我需要定义与该像素对应的区域(25 * 25 公里或 0.25 度)。为此，我必须知道该像素(单元格)的 4 个坐标(lat-long)。 感谢您的帮助

Answer 1

EASE GRID 使用公制 EPSG 3410 定义的全局圆柱等积投影。只要我看到它，空间范围应该以米为单位提供，而不是地理坐标。 来自 here我们看到 map 范围坐标是:

最小值:-17609785.303313

ymin: -7389030.516717

xmax: 17698276.686747

ymax: 7300539.133283

稍微改变一下你的代码，我们就有了这个

library(raster)
library('rgdal')
wdata <- 'D:/Programacao/R/Raster/Stackoverflow'
wshp <- 'S:/Vetor/Administrativo/Portugal'
#setwd(wdata)
file <- readBin(file.path(wdata, "ID2r1-AMSRE-ML2010001D.v03.06H"),
               integer(), size=2,  n=586 * 1383, signed=T)
m <- matrix(data = file, ncol = 1383, nrow = 586, byrow = TRUE)
-17609785.303313 -7389030.516717 17698276.686747 7300539.133283
rm <- raster(m, xmn = -17609785.303313, xmx = 17698276.686747,
             ymn = -7389030.516717, ymx = 7300539.133283)
proj4string(rm) <- CRS('+init=epsg:3410')

> rm
class       : RasterLayer 
dimensions  : 586, 1383, 810438  (nrow, ncol, ncell)
resolution  : 25530.05, 25067.53  (x, y)
extent      : -17609785, 17698277, -7389031, 7300539  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3410 +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs 
data source : in memory
names       : layer 
values      : 0, 3194  (min, max)

writeRaster(rm, file.path('S:/Temporarios', 'easegrridtest.tif'), overwrite = TRUE)
plot(rm, asp = 1)

我们现在可以叠加一些空间数据

afr <- readOGR(dsn = file.path(wdata), layer = 'Africa_final1_dd84')
proj4string(afr) <- CRS('+init=epsg:4326') # Asign projection 
afr1 <- spTransform(afr, CRS(proj4string(rm)))

plot(afr1, add = T)

现在您可以开始使用 ROI 的提取范围，可能使用 extent()

我对空间调整不满意。通过这种方法，我遇到了一个巨大的错误。我不确定位置错误，但针对此产品进行了描述。也许有范围参数。

由于您有兴趣将其用于地面测量，您可以将您的 ROI 多边形化并在您的 GPS 或 GIS 中使用它。

您还可以通过以下方式获得感兴趣的细胞范围:

选择一个近似坐标，识别像元并获取范围:

cell <- cellFromXY(rm, matrix(c('x'= -150000, 'y' =200000), nrow = 1, byrow = T))
r2 <- rasterFromCells(rm, cell, values=TRUE)
extent(r2)
class       : Extent 
xmin        : -172759.8 
xmax        : -147229.7 
ymin        : 181362 
ymax        : 206429.6 

并可能在 map (或绘图)中识别 ROI(单个单元格)

cell <- cellFromXY(rm, matrix(c('x'= -1538000, 'y' =1748000), nrow = 1, byrow = T))
r2 <- rasterFromCells(rm, cell, values=TRUE)
r2p <- as(r2, 'SpatialPolygons')
extr2 <- extent(r2) + 300000
plot(rm,  col = heat.colors(6), axes = T, ext = extr2)
plot(afr1, add = T, col = 'grey70')
plot(r2p, add = T)

针对特定区域

假设单元格是指列值和行值，您可以继续 raster::cellFromRowCol

cell2 <- cellFromRowCol(rm, rownr = mycellnrow, colnr = mycellnrow)
r3 <- rasterFromCells(rm, cell2, values=TRUE)
r3p <- as(r3, 'SpatialPolygons')
extr3 <- extent(r3) + 3000000

在这种特殊情况下，123 和 450 似乎远离任何大陆区域......

希望对您有所帮助。

有关 AMSR-E/Aqua 每日网格亮度温度的更多信息 here