我如何正确地从 R 中的 NETCDF 文件中重新定位 RasterBrick 中错位的经度和纬度？
How do I, correctly, reposition misplaced longitude and latitude in a RasterBrick from NETCDF file in R?
我有一个 RasterBrick 包含来自 TRMM (TMPA/3B43) Rainfall Estimate L3 1 month 0.25 degree x 0.25 degree V7 (TRMM_3B43. Extract from data is in the tif file 的降水量变量。
此 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 。
图片实际应该如 from Image source
所示
我尝试使用 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。所以，如果你也知道的话，我将不胜感激。
谢谢
这是一个简单的转置和翻转数据的工作流程；并设置正确的范围。我从您指向的源下载了一个文件，但这些是 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)
我的理解是第一层是“沉淀”（mm/hr），第二层是“相对误差”，第三层是“gaugeRelativeWeighting”。
如果你必须多次这样做，你可以使用这样的函数（它也可以抑制警告；并且只 returns 第一层）；并使用 |> 你需要 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是lon/lat（这不是MODIS数据）
这就是我能够解决数据的单变量版本的方法。但是，在此过程中删除了日期戳。
#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)
我有一个 RasterBrick 包含来自 TRMM (TMPA/3B43) Rainfall Estimate L3 1 month 0.25 degree x 0.25 degree V7 (TRMM_3B43. Extract from data is in the tif file 的降水量变量。
此 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。所以，如果你也知道的话，我将不胜感激。
谢谢
这是一个简单的转置和翻转数据的工作流程；并设置正确的范围。我从您指向的源下载了一个文件，但这些是 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)
我的理解是第一层是“沉淀”（mm/hr），第二层是“相对误差”，第三层是“gaugeRelativeWeighting”。
如果你必须多次这样做，你可以使用这样的函数（它也可以抑制警告；并且只 returns 第一层）；并使用 |> 你需要 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是lon/lat（这不是MODIS数据）
这就是我能够解决数据的单变量版本的方法。但是，在此过程中删除了日期戳。
#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)