将 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)

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

Converting HDF to georeferenced file (geotiff, shapefile)

r

raster

hdf