如何正确地从R中的RasterBrick文件中重新定位经度和纬度？

Question

我有一个RasterBrick，包含来自3B43的降水变量。从数据中提取的内容在tif文件中。

此RasterBrick是使用以下代码从NETCDF file生成的

library(raster)

# set working directory    
setwd("C:/_data/TRMM_data/TRMM")

# import netcdf file into R using brick function in raster package
ppt.brick <- brick("TRMM_3B43_Aggregation.ncml.ncml.nc4")

# set coordinate resference system
crs(ppt.brick) = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m"

数据ppt.brick提供了如下信息：

> ppt.brick
class      : RasterBrick 
dimensions : 1440, 400, 576000, 252  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -50, 50, -180, 180  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
source     : C:/sijeh_data/TRMM_data/TRMM/TRMM_3B43_Aggregation.ncml.ncml.nc4 
names      : X1999.01.16, X1999.02.15, X1999.03.16, X1999.04.16, X1999.05.16, X1999.06.16, X1999.07.16, X1999.08.16, X1999.09.16, X1999.10.16, X1999.11.16, X1999.12.16, X2000.01.16, X2000.02.15, X2000.03.16, ... 
Date       : 1999-01-16, 2019-12-16 (min, max)
varname    : precipitation

plot(ppt.brick)生成下面的image。但是正如你所看到的，latitude在x-axis而不是y-axis上绘制，longitude在y-axis上而不是x-axis上。

。

该图像实际上应该如

来自图像源

我尝试过使用transpose函数ppt.t <- t(ppt.brick)，但这会导致RasterBrick丢失一些信息，比如varname、date和time，如下所示。

>ppt.t
class      : RasterBrick 
dimensions : 400, 1440, 576000, 252  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs 
source     : C:/Users/saa815/AppData/Local/Temp/RtmpIHEl9R/raster/r_tmp_2021-05-19_141009_7792_98365.grd 
names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,  layer.7,  layer.8,  layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
min values :        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0, ... 
max values : 1.837709, 1.792068, 2.171153, 1.309954, 1.687151, 2.159085, 3.079944, 2.278921, 1.547583, 1.631066, 2.687030, 1.929542, 1.488750, 1.950469, 1.307860, ...

请问什么是最好的方法，以transpose一个RasterBrick，但仍保留在原始data的信息？也很难判断transpose是否正确，因为应该是-50的轴在transpose之后变成+50。见图

​

之后，我还需要使用一个mask RasterBrick和一个polygon。所以，如果你也知道的话，我会很感激的。

谢谢

Answer 1

下面是一个简单的工作流，它可以转换和翻转数据，并设置正确的范围。我从你指向的源代码下载了一个文件，但这些文件是HDF文件，而不是ncdf文件。否则一切似乎都是一样的。

library(terra)
#terra version 1.2.13

f <- "3B43.20160701.7.HDF"
x <- rast(f)
#Warning message:
#[rast] unknown extent

y <- t(x)
z <- flip(y, "v")
ext(z) <- c(-180, 180, -50, 50)

z
#class       : SpatRaster 
#dimensions  : 400, 1440, 3  (nrow, ncol, nlyr)
#resolution  : 0.25, 0.25  (x, y)
#extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#names       :           0,           1,           2 
#min values  : 0.000000000, 0.004333077, 0.000000000 
#max values  :   2.0519991,   0.5469698,  98.0000000 

plot(z, nc=1)

​

​

我的理解是，第一层是“降水量”(毫米/小时)，第二层是“相对误差”，第三层是"gaugeRelativeWeighting“。

如果您必须多次这样做，您可以使用类似这样的函数(这也会抑制警告，并且只返回第一层)；以及使用需要R版本4.1的|>。

get_3B43 <- function(filename) {
   on.exit(options(warn=options()$warn))
   options(warn=-1)
   subf <- paste0('HDF4_SDS:UNKNOWN:"', filename, '":0"')
   x <- rast(subf) |> t() |> flip(direction="v")
   ext(x) <- c(-180, 180, -50, 50)
   time(x) <- as.Date(substr(filename, 6, 13), "%Y%m%d")
   names(x) <- "precipitation"
   units(x) <- "mm/hr" 
   x
}

r <- get_3B43(f)

r
#class       : SpatRaster 
#dimensions  : 400, 1440, 1  (nrow, ncol, nlyr)
#resolution  : 0.25, 0.25  (x, y)
#extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
#source      : memory 
#name        : precipitation 
#min value   :             0 
#max value   :      2.051999 
#unit        :         mm/hr 
#time        : 2016-07-01 

有了你提供的文件，我可以

ff <- list.files()
x <- lapply(ff, get_3B43)
x <- rast(x)
x
class       : SpatRaster 
dimensions  : 400, 1440, 12  (nrow, ncol, nlyr)
resolution  : 0.25, 0.25  (x, y)
extent      : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
sources     : memory  
              memory  
              memory  
              ... and 9 more source(s)
names       : preci~ation, preci~ation, preci~ation, preci~ation, preci~ation, preci~ation, ... 
min values  :           0,           0,           0,           0,           0,           0, ... 
max values  :    1.837709,    1.792068,    2.171153,    1.309954,    1.687151,    2.159085, ... 
unit        :       mm/hr,       mm/hr,       mm/hr,       mm/hr,       mm/hr,       mm/hr, ... 
time        : 1999-01-01 to 1999-12-01 

您可能需要成批地这样做(例如，一年一次，然后将它保存到磁盘上)，也许是这样

z <- writeCDF(x, "test.nc", varname="precipitation", unit="mm/hr")

另外，我不知道你为什么这样做：

crs(ppt.brick) = "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m"

这是错误的，CRS是龙/拉特(这不是MODIS的数据)。

Answer 2

这就是我如何能够解决单一变量版本的数据。但是，在此过程中删除了日期标记。

#Required package 

library(raster)

#list nc files in the path
lf = list.files('your path', pattern = ".nc$")

#stack the files
s = stack(lf)

#convert to RasterBrick
b = brick(s)

#assign CRS
crs(b) = "+proj=longlat +datum=WGS84 +no_defs"

#confirm the extent
extent(b) = c(-50, 50,-180, 180)

#flip transpose and flip
bb = flip(t(flip(b,1)), 1)

得到的RasterBrick如下所示

> bb
class      : RasterBrick 
dimensions : 400, 1440, 576000, 252  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : -180, 180, -50, 50  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : C:/Users//AppData/Local/Temp/RtmpgZPSMv/raster/r_tmp_2021-05-20_163037_14372_39367.grd 
names      :  layer.1,  layer.2,  layer.3,  layer.4,  layer.5,  layer.6,  layer.7,  layer.8,  layer.9, layer.10, layer.11, layer.12, layer.13, layer.14, layer.15, ... 
min values :        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0,        0, ... 
max values : 1.837709, 1.792068, 2.171153, 1.309954, 1.687151, 2.159085, 3.079944, 2.278921, 1.547583, 1.631066, 2.687030, 1.929542, 1.488750, 1.950469, 1.307860, ... 

图像正确地对齐，如图所示

plot(bb[[1]])

​

​

或者，如果使用NetCDF或NetCDF4文件，也可以这样做。

library(ncdf4)
library(raster)

f = "TRMM_3B43_Aggregation.ncml.ncml.nc4"

ncs <- nc_open(f)
print(ncs)
precip = ncvar_get(ncs, "precipitation")
r=raster(precip[,,1]) 
crs(r) = "+proj=longlat +datum=WGS84 +no_defs"
extent(r) = c(-180, 180,-50,50)
r = flip(r,c(2))
plot(r)