将HDF转换为地理参考文件(geotiff，shapefile) [英] Converting HDF to georeferenced file (geotiff, shapefile)

问题描述

我正在处理28个关于海洋主要生产力的HDF4文件(年度.tar文件可在此处找到:

I am dealing with 28 HDF4 files on ocean primary productivity (annual .tar files can be found here: http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.cbpm2.v.php) My goal is to do some calculations (I need to calculate concentrations per area and obtain the mean over several years, i.e. combine all files spatially) and then convert them to a georeferenced file I can work with in ArcGIS (preferably shapefile, or geotiff).

我尝试了几种方法来转换为ASCII或光栅文件，然后使用gdalUtils工具(例如gdal_translate和get_subdatasets)添加投影.但是，由于HDF4文件没有按标准命名(与MODIS文件不同)，因此后者不起作用，并且我无法访问这些子集.

I have tried several ways to convert to ASCII or raster files and then add a projection using gdalUtils tools such as gdal_translate and get_subdatasets. However, as the HDF4 files are not named after the standards (unlike the MODIS files), the latter doesn't work and I cannot access the subsets.

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

Here's the code I used to convert to raster:

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)

生成的投影完全关闭.我非常感谢您提供帮助(最好通过批处理的方式来做到这一点(通过任何有效的格式进行转换).

The resulting projection is completely off. I'd appreciate help how to do this (converting through any format that works), preferably as batch process.

推荐答案

您尚未设置栅格的范围，因此假设栅格的比例为1:ncols，1:nrows，这对经纬度的定位不正确数据集...

You've not set the extent of your raster, so its assumed to be 1:ncols, 1:nrows, and that's not right for a lat-long data set...

gdalinfo表示其意图是一个完整的球体，所以如果我这样做:

The gdalinfo implies its meant to be a full sphere, so if I do:

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

我看到的栅格具有完整的全球纬度-纬度范围，当我将栅格加载到QGIS中时，它会与OpenStreetMap很好地叠加.

I see a raster with full global lat-long extent and when I load the raster into QGIS it overlays nicely with an OpenStreetMap.

文件中似乎没有足够的元数据来进行精确的投影(地球半径和偏心率是多少?)，所以不要尝试对此数据进行任何小规模的操作...

There doesn't seem to be enough metadata in the file to do a precise projection (what's the Earth radius and eccentricity?) so don't try and do anything small-scale with this data...

这是它的外观:

此外，您还跳过了一些不必要的阅读框.您可以直接阅读HDF并设置其范围和投影:

Also you've jumped through a few unnecessary hoops to read this. You can read the HDF directly and set its extent and projection:

> 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:

And land values to NA:

> r[r==-9999]=NA

写出来，画出来:

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

这篇关于将HDF转换为地理参考文件(geotiff，shapefile)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！