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 栅格像元的坐标。
我有一个使用 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 栅格像元的坐标。