R如何使用砖块从常规tif中获取纬度向量

R how to get a vector of latitudes from a regular tif using brick

我有一个使用 brick 读取的 tif 文件。我可以看到 tif 文件的空间扩展:

b<-brick("t.tif")
b 
class       : RasterBrick 
dimensions  : 10, 10, 100, 1  (nrow, ncol, ncell, nlayers)
resolution  : 0.0001851853, 0.0001851854  (x, y)
extent      : -18.61944, -18.61759, 37.83856, 37.84041  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /scratch/tompkins/water_allafrica/t.tif 
names       :   t 
min values  :   0 
max values  : 255 

extent(SpatialPoints(b))
class       : Extent 
xmin        : -18.61935 
xmax        : -18.61769 
ymin        : 37.83865 
ymax        : 37.84031 

但我无法弄清楚如何轻松获得纬度和经度矢量,我需要定义一个我想写出来的 netcdf 文件头。我可以手动完成,但我假设有一个更易于使用的内置函数。

示例输入文件在这里:http://clima-dods.ictp.it/Users/tompkins/Whosebug/t.tif

这可能是您需要的:

ext <- extent(b)

lat <- seq(ext@ymin,ext@ymax,res(b)[2])
lon <- seq(ext@xmin,ext@xmax,res(b)[1])

所以基本上您是在创建一个从 x/y 最小值到最大值的序列向量,其中包含积木分辨率的间距。

这些值是指单元格的角坐标...您可能也对单元格中心感兴趣。

仅供说明:

# create testraster
x <- raster(resolution=c(40,40))

x[]<- 1:ncell(x)

# plot

plot(x)

# add corner coordinates


plot(SpatialPoints(cbind(rep(extent(x)@xmin,10),seq(extent(x)@ymin,extent(x)@ymax,res(x)[2])),proj4string = crs(x)),
     col='red',pch='*',cex=5,add=T)


# add cell centers

plot(SpatialPoints(xyFromCell(x,cellFromRowCol(x,1:nrow(x),1)),proj4string = crs(x)),
     col='blue',pch='*',cex=5,add=T)

所以上面的方法给出了红色星号指示的纬度。如果你需要蓝色的,你可以使用 xyFromCell 其中 returns 栅格像元的坐标。