Преобразование HDF в файл с географической привязкой (geotiff, shapefile)

Question

Я имею дело с 28 файлами HDF4 по первичной продуктивности океана (годовые файлы.tar можно найти здесь: http://orca.science.oregonstate.edu/1080.by.2160.monthly.hdf.cbpm2.v.php). Цель состоит в том, чтобы сделать некоторые расчеты (мне нужно рассчитать концентрации на площадь и получить среднее значение за несколько лет, т.е. объединить все файлы пространственно), а затем преобразовать их в файл с географической привязкой, с которым я могу работать в ArcGIS (предпочтительно, шейп-файл или геотиф).

Я пробовал несколько способов конвертировать в 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, и это неправильно для набора данных длиной в lat...

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"

И стоимость земли для НС:

> r[r==-9999]=NA

Напишите это, постройте это:

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