[R-sig-Geo] Converting .hdf file to raster layer
Michael Sumner
mdsumner at gmail.com
Thu Aug 14 10:52:30 CEST 2014
Just had a look, the file has subdatasets and seemingly raster can't
handle those. (I haven't explored this much but might have a chance
to).
library(raster)
library(rgdal) ## built with HDF4, HDF5 and NetCDF4
fsrc <- "ftp://ftp.star.nesdis.noaa.gov/pub/corp/scsb/wguo/data/VHP_16km/VH/VHP.G16.C07.NC.P1981035.ND.hdf"
f <- basename(fsrc)
if (!file.exists(f)) download.file(fsrc, f, mode = "wb")
## found with system GDAL, see below
sds <- c('HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":0',
'HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":1',
'HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":2'
)
##r <- brick(sds)
##Error in GDALinfo(filename, silent = silent, returnRAT = RAT,
returnCategoryNames = RAT) :
## object 'RATlist' not found
##Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer", :
## Cannot create a RasterLayer object
from this file.
## In addition: Warning message:
## In dim(x) : no bands in dataset
## try with rgdal directly, so far so good but
## no spatial-reference (again see gdalinfo output below)
raster(readGDAL(sds[1]) )
raster(readGDAL(sds[1]) )
HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":0 has GDAL driver HDF4Image
and has 904 rows and 2500 columns
class : RasterLayer
dimensions : 904, 2500, 2260000 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 0, 2500, 0, 904 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : in memory
names : band1
values : -99.9, 322 (min, max)
Warning message:
In readGDAL(sds[1]) : GeoTransform values not available
## proceed without caution, again values from gdalinfo
b <- brick(stack(sapply(sds, function(x) raster(readGDAL(x)))))
proj4string(b) <- CRS(" +proj=longlat +ellps=WGS84 +datum=WGS84
+no_defs +towgs84=0,0,0")
extent(b) <- c(-180, 180, -55.15200043, 75.02400208)
library(maptools
data(wrld_simpl)
plot(b, addfun = function() plot(wrld_simpl, add = TRUE))
That looks ok to me, I don't know if it's WGS84 or the sphere or
something else.
Cheers, Mike.
system(sprintf("gdalinfo %s", f))
Driver: HDF4/Hierarchical Data Format Release 4
Files: VHP.G16.C07.NC.P1981035.ND.hdf
Size is 512, 512
Coordinate System is `'
Metadata:
ANCILLARY_FILES=FILE_CONFIGURE:vh.config
FILE_PRELAUNCH_CALIBRATION:../ancillary/AVHRR_calibration_prelaunch.txt
FILE_POSTLAUNCH_CALIBRATION:../ancillary/AVHRR_calibration_postlaunch.txt
FILE_IGBP_LANDTYPE:../ancillary/igbp_landtype_GVIx.hdf
FILE_METADATA_REGIONS:../ancillary/regions_for_metadata.txt
FILE_EDF_NDVI:../ancillary/NVI_counts_ByLine_G04.hdf
FILE_EDF_BT4:../ancillary/BT4_counts_ByLine_G04.hdf
CITATION_TO_DOCUMENTS=User Guide of Vegetation Health(VH) system
(version 1.3, March 21 2012)
CONFIGURE_FILE_CONTENT=[Options for vh.exe]
DIR_Ancillary= ../ancillary
DIR_GVI= data/VH_unitTest1/weekly
DIR_VH= data/VH_unitTest1/G04
DIR_VH_META= data/VH_unitTest1/G04/meta
DIR_CLIMAT= data/VH_EDF_v2/climate_G04
DIR_CLIMAT_META= data/VH_EDF_v2/climate_G04/meta
FILE_PREFIX= VHP
Input_Data_Type= VHP
ResolutionString= G04
Days_Per_Period= 7
FilterSize= 15
applyEDFonNDVI= 1 # 0: none, 1: linebyline for NVI;
applyEDFonBT= 1
Instrument= AVHRR # so far AVHRR is the only option
FORMAT_GVI= NETCDF #or HDF4
FORMAT_CLIMAT= HDF4
FORMAT_ND= NETCDF
FORMAT_SM= NETCDF
FORMAT_VH= NETCDF
[Periods of GVI data used for VH]
# this section control which satellite will be used for calculating
ND, SM and VH
#satID satNumber yearWeek1 yearWeek2
NC 07 198135 198449
NF 09 198509 198844
NH 11 198846 199436
NJ 14 199504 200052
NL 16 200101 200401
NL 16 200405 200410
NL 16 200425 200428
NL 16 200430 200523
NN 18 200524 201052
NP 19 201101 399999
[Periods of AVHRR data used for GVI climatology]
#this section controls which satellite will be used for creating VH climatology
#satID satNumber yearWeek1 yearWeek2
NC 07 198142 198450
NF 09 198515 198752
NH 11 198920 199252
NJ 14 199520 199952
NL 16 200120 200252
[END]
CONTACT=NOAA/NESDIS/STAR/EMB
DATE_BEGIN=239
DATE_END=245
DAYS_PER_PERIOD=7
END_LATITUDE_RANGE=-55.15200043
END_LONGITUDE_RANGE=180
INPUT_FILENAMES=data/VH_unitTest1/weekly/VHP.G04.C07.NC.P1981035.S239.E245.nc
INPUT_FILES=1
INSTRUMENT=AVHRR
PERIOD_OF_YEAR=35
PRODUCT_NAME=Vegetation Health
PROJECTION=Plate_Carree
SATELLITE=NC
START_LATITUDE_RANGE=75.02400208
START_LONGITUDE_RANGE=-180
TIME_BEGIN=00:00 UTC (use day time data only)
TIME_END=23:59 UTC (use day time data only)
VERSION=VH (vh.exe,version 1.3, March 21 2012)
YEAR=1981
Subdatasets:
SUBDATASET_1_NAME=HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":0
SUBDATASET_1_DESC=[904x2500] BT4 (16-bit integer)
SUBDATASET_2_NAME=HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":1
SUBDATASET_2_DESC=[904x2500] NDVI (16-bit integer)
SUBDATASET_3_NAME=HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":2
SUBDATASET_3_DESC=[904x2500] QA (8-bit unsigned integer)
Corner Coordinates:
Upper Left ( 0.0, 0.0)
Lower Left ( 0.0, 512.0)
Upper Right ( 512.0, 0.0)
Lower Right ( 512.0, 512.0)
Center ( 256.0, 256.0)
On Thu, Aug 14, 2014 at 10:27 AM, Arnold Salvacion
<arnold_salvacion at yahoo.com> wrote:
> Dear Colleagues,
>
> Does any here have already experience converting NOAA AVHRR VHP Product in .hdf (ftp://ftp.star.nesdis.noaa.gov/pub/corp/scsb/wguo/data/VHP_16km/VH/ ) format to a raster layer in R?
>
>
> Best regards,
>
> Arnold
> [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Michael Sumner
Software and Database Engineer
Australian Antarctic Division
Hobart, Australia
e-mail: mdsumner at gmail.com
More information about the R-sig-Geo
mailing list