Преобразование HDF4 в растр; с сетками долготы/широты доступны в другом файле HDF

Я пытаюсь преобразовать файлы HDF4 (представляющие ежедневную концентрацию морского льда) в растровый объект в R. Однако сами файлы HDF не содержат сетки долготы/широты или информацию о проекции, и такая информация должна быть извлечена из другого файла hdf.

На веб-сайте о формате данных говорится:

Data Format
Sea ice concentration maps with two different color scales are available as PNG image. The NIC color scale uses the same colors as the National Ice Center, the "visual" color scale uses white and shades of grey.
There is one file per day per region per color scale.
Sea ice concentration data are available as HDF4 files: There is one file per day per region. Each file contains one two-dimensional array of the sea ice concentration in a polar stereographic grid.
The longitude and latitude coordinates of each pixel in a the HDF4 file are saved in extra files, one file per region for each available resolution.
They are found here: https://seaice.uni-bremen.de/data/grid_coordinates/,  sorted by hemisphere and grid resolution (see also the README file https://seaice.uni-bremen.de/data/grid_coordinates/README).
GEOTIFF files use the NIC color scale and were tested to work with QGIS. Ice concentrations are scaled between 0 and 100, land and missing values are set to 120 (older files: SIC: 0-200, land/NaN: 255).

Я пытался использовать R для загрузки эта карта, используя этот код:

> require(raster)
> CurrTemp <- tempfile()
> download.file(url = "https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf", destfile = CurrTemp, mode = "wb", quiet = T)
> Map1 <- readAll(raster(CurrTemp))
> plot(Map1)
> Map1
class       : RasterLayer 
dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 1264, 0, 1328  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : in memory
names       : file43fc5b4e68de 
values      : 0, 100  (min, max)

Карта загружается в R как растровый объект, но с неправильными координатами и без проекции. Согласно этой странице, координаты должны быть извлечены из еще один файл hdf.

Не могли бы вы сообщить мне, как преобразовать эти файлы hdf в растровые объекты с правильными координатами и проекцией.

Спасибо.


r gis rgdal r-raster hdf
person Ahmed El-Gabbas    schedule 21.09.2018    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Я использовал один из файлов geotiff, которые они также предоставляют, чтобы найти экстент и crs.

library(raster)
raster('asi-AMSR2-s6250-20180922-v5.tif')
#class       : RasterLayer 
#dimensions  : 1328, 1264, 1678592  (nrow, ncol, ncell)
#resolution  : 6250, 6250  (x, y)
#extent      : -3950000, 3950000, -3950000, 4350000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs 
#data source : asi-AMSR2-s6250-20180922-v5.tif 
#names       : asi.AMSR2.s6250.20180922.v5 
#values      : 0, 255  (min, max)

Теперь я знаю, что могу сделать

library(raster)
CurrTemp <- tempfile()
download.file(url = "https://seaice.uni-bremen.de/data/amsre/asi_daygrid_swath/s6250/2003/feb/Antarctic/asi-s6250-20030214-v5.hdf", destfile = CurrTemp, mode = "wb", quiet = T)
r <- raster(CurrTemp)

extent(r) <- c(-3950000, 3950000, -3950000, 4350000)
crs(r) <- "+proj=stere +lat_0=-90 +lat_ts=-70 +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs "
# writeRaster(r, 'my_asi-s6250-20030214-v5.tif')

Файл «другой hdf» имеет значения долготы/широты для ячеек, но это не то, что вам нужно, поскольку данные не имеют системы координат долготы/широты.

person Robert Hijmans    schedule 23.09.2018