将 HDF4 转换为光栅; longitude/latitude 网格在另一个 HDF 文件中可用

Convert HDF4 to raster; with longitude/latitude grids are available in another HDF file

我正在尝试将 HDF4 文件(代表每日海冰浓度)转换为 R 中的栅格对象。但是,HDF 文件本身不包含 longitude/latitude 网格或投影信息,此类信息应该是从另一个 hdf 文件中提取。

关于数据格式的website说:

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 加载 this map 使用此代码:

> 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 中,但坐标错误且没有投影。根据this page, coordinates should be extracted from another hdf file.

能否请您告诉我如何将这些 hdf 文件转换为具有正确坐标和投影的光栅对象。

谢谢。

我使用了他们提供的其中一个 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')

"other hdf" 文件具有单元格的经度/纬度值,但这不是您想要的,因为数据没有 lon/lat 坐标参考系统。