将 HDF 转换为地理参考文件（geotiff、shapefile）

Question

我正在处理关于海洋初级生产力的 28 个 HDF4 文件（可以在此处找到年度 .tar 文件：http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.cbpm2.v.php） 我的目标是进行一些计算（我需要计算每个区域的浓度并获得几年的平均值，即在空间上合并所有文件），然后将它们转换为我可以在 ArcGIS 中使用的地理参考文件（最好是 shapefile 或 geotiff） .

我尝试了几种方法来转换为 ASCII 或光栅文件，然后使用 gdalUtils 工具（例如 gdal_translate 和 get_subdatasets）添加投影。但是，由于 HDF4 文件不是以标准命名的（与 MODIS 文件不同），因此后者不起作用，我无法访问子集。

这是我用来转换为栅格的代码：

library(raster)
library(gdalUtils)

setwd("...path_to_files...")

gdalinfo("cbpm.2015060.hdf")
hdf_file <- "cbpm.2015060.hdf"

outfile="testout"
gdal_translate(hdf_file,outfile,sds=TRUE,verbose=TRUE)
file.rename(outfile,paste("CBPM_test",".tif",sep="")) 

rast <- raster("CBPM_test.tif")

wgs1984 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
projection(rast) <- wgs1984
#crs(rast) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 

plot(rast)

writeRaster(rast, file="CBPM_geo.tif", format='GTiff', overwrite=TRUE)

生成的投影完全关闭。我会很感激帮助如何做到这一点（通过任何有效的格式转换），最好是批处理。

Answer 1

您尚未设置栅格的范围，因此假定其为 1:ncols、1:nrows，这不适用于经纬度数据集...

gdalinfo 暗示它是一个完整的球体，所以如果我这样做：

 extent(rast)=c(xmn=-180, xmx=180, ymn=-90, ymx=90)
 plot(rast)
 writeRaster(rast, "output.tif")

我看到一个具有完整全球经纬度范围的栅格，当我将栅格加载到 QGIS 中时，它与 OpenStreetMap 很好地重叠。

文件中似乎没有足够的元数据来进行精确投影（地球半径和偏心率是多少？）所以不要尝试用这些数据做任何小规模的事情......

它的外观如下：

您还跳过了一些不必要的障碍来阅读本文。您可以直接读取 HDF 并设置其范围和投影：

> r = raster("./cbpm.2017001.hdf")

我们得到了什么：

> r
class       : RasterLayer 
dimensions  : 1080, 2160, 2332800  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 2160, 0, 1080  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : /home/rowlings/Downloads/HDF/cbpm.2017001.hdf 
names       : cbpm.2017001 

设置范围：

> extent(r)=c(xmn=-180, xmx=180, ymn=-90, ymx=90)

和投影：

> projection(r)="+init=epsg:4326"

以及对 NA 的土地价值：

> r[r==-9999]=NA

写下来，画出来：

> writeRaster(r,"r.tif")
> plot(r)