Конвертировать HDF4 в растр;с сетками долготы / широты доступны в другом файле HDF - PullRequest
Новая
       105

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

0 голосов
/ 21 сентября 2018

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

Веб-сайт в формате данных сообщает: 

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-файлы в растровые объекты с правильными координатамии проекция.

Спасибо.

1 Ответ

0 голосов
/ 23 сентября 2018

Я использовал один из файлов 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""В файле есть значения долготы / широты для ячеек, но это не то, что вам нужно, поскольку у данных нет системы координат lon / lat.

...