代码之家 › 专栏 › 技术社区 › Ahmed El-Gabbas

将HDF4转换为光栅;在另一个HDF文件中提供了经度/纬度网格

hdf rgdal r-raster gis r

1

Ahmed El-Gabbas · 技术社区 · 7 年前

我正在尝试将HDF4文件(表示每日海冰浓度)转换为R中的光栅对象。但是,HDF文件本身不包含经纬度网格或投影信息,这些信息应该从另一个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 another hdf file .

1 回复 | 直到 7 年前

1

Robert Hijmans 7 年前

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坐标系。