[R-sig-Geo] Converting .hdf file to raster layer

Arnold Salvacion arnold_salvacion at yahoo.com
Thu Aug 14 17:30:19 CEST 2014 

Hi Mike,

Thanks. I tried you code  but got the error message below:

Arnold


R version 3.0.3 (2014-03-06) -- "Warm Puppy"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R. > 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") trying URL 'ftp://ftp.star.nesdis.noaa.gov/pub/corp/scsb/wguo/data/VHP_16km/VH/VHP.G16.C07.NC.P1981035.ND.hdf' using Synchronous WinInet calls opened URL downloaded 2.1 Mb > 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' + ) > raster(readGDAL(sds[1])  ) Error in raster(readGDAL(sds[1])) :  error in evaluating the argument 'x' in selecting a method for function 'raster': Error in .local(.Object, ...) :  `C:\Users\Owner\Documents\HDF4_SDS:UNKNOWN:"VHP.G16.C07.NC.P1981035.ND.hdf":0' does not exist in the file system,
and is not recognised as a supported dataset name. 



On Thursday, August 14, 2014 4:52 PM, Michael Sumner <mdsumner at gmail.com> wrote:
 


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
	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list