[R-sig-Geo] Converting .hdf file to raster layer
Ahmadou Dicko
dicko.ahmadou at gmail.com
Thu Aug 21 11:31:34 CEST 2014
I suspect that your GDAL was not build with the HDF4 driver.
What is the output of :
grep("hdf4", gdalDrivers()$name, ignore.case = TRUE, value = TRUE)
On Thu, Aug 21, 2014 at 7:52 AM, Arnold Salvacion <
arnold_salvacion at yahoo.com> wrote:
> Hi Ahmadou,
>
> Have tried your code but still no success. However, I have actually tried
> to download the .hdf file directly from net. What I did is import the data
> from local folder. Below is the error that I got:
>
> Error in .local(.Object, ...) :
> `C:\Users\Owner\Documents\VHP.G16.C07.NC.P1982006.VH.hdf' not recognised as a supported file format.
>
>
> Best regards,
>
> Arnold
>
>
>
>
> On Friday, August 15, 2014 2:33 AM, Ahmadou Dicko <
> dicko.ahmadou at gmail.com> wrote:
>
>
> You can use the gdalUtils package to get information on datasets inside
> your hdf file and then stack them using the raster package
>
> library(raster)
> library(rgdal) ## built with HDF4 support
> library(gdalUtils)
>
> fsrc <- "
> ftp://ftp.star.nesdis.noaa.gov/pub/corp/scsb/wguo/data/VHP_16km/VH/VHP.G16.C07.NC.P1981035.ND.hdf
> "
> f <- file.path("/tmp", basename(fsrc))
> if (!file.exists(f)) download.file(fsrc, f, mode = "wb")
>
> info <- gdalinfo(f)
> data <- grep("HDF4_SDS", info, value = TRUE)
> data <- gsub("SUBDATASET_\\d+_NAME=|\\s+", "", data)
>
>
> sds <- stack(data)
> sds
> ## class : RasterStack
> ## dimensions : 904, 2500, 2260000, 3 (nrow, ncol, ncell, nlayers)
> ## resolution : 1, 1 (x, y)
> ## extent : 0, 2500, 0, 904 (xmin, xmax, ymin, ymax)
> ## coord. ref. : NA
> ## names : VHP.G16.C07.NC.P1981035.ND.1,
> VHP.G16.C07.NC.P1981035.ND.2, VHP.G16.C07.NC.P1981035.ND.3
> ## min values : -32768,
> -32768, 0
> ## max values : 32767,
> 32767, 255
>
> ##You can use RasterBrick too, by doing `brick(sds)`
>
> sessionInfo()
> ## R version 3.1.1 Patched (2014-08-12 r66349)
> ## Platform: x86_64-unknown-linux-gnu (64-bit)
>
> ## locale:
> ## [1] LC_CTYPE=en_US.utf8 LC_NUMERIC=C
> ## [3] LC_TIME=en_US.utf8 LC_COLLATE=en_US.utf8
> ## [5] LC_MONETARY=en_US.utf8 LC_MESSAGES=en_US.utf8
> ## [7] LC_PAPER=en_US.utf8 LC_NAME=C
> ## [9] LC_ADDRESS=C LC_TELEPHONE=C
> ## [11] LC_MEASUREMENT=en_US.utf8 LC_IDENTIFICATION=C
>
> ## attached base packages:
> ## [1] stats graphics grDevices utils datasets methods
> ## [7] base
>
> ## other attached packages:
> ## [1] gdalUtils_0.3.2 rgdal_0.8-16 raster_2.2-32
> ## [4] sp_1.0-15
>
> ## loaded via a namespace (and not attached):
> ## [1] codetools_0.2-8 compiler_3.1.1 foreach_1.4.2
> ## [4] grid_3.1.1 iterators_1.0.7 lattice_0.20-29
> ## [7] R.methodsS3_1.6.1 R.oo_1.18.0 R.utils_1.32.4
> ## [10] tools_3.1.1
>
> grep("hdf4", gdalDrivers()$name, ignore.case = TRUE, value = TRUE)
> ## [1] "HDF4" "HDF4Image"
>
>
>
> On Thu, Aug 14, 2014 at 8:52 AM, 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
> <http://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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
>
>
> --
> Ahmadou H. DICKO
> statistician economist (Ingénieur Statisticien Économiste)
> PhD candidate in Climate change economics
> Faculty of economics and managment - Cheikh Anta Diop University
> West African Science Service Center on Climate Change and Adaptated Land
> Use (WASCAL)
> Center for Development Research (ZEF) - University of Bonn
>
> email :
> ahmadou.dicko at ucad.edu.sn
> twitter : @dickoah
> github : github/dickoa <https://github.com/dickoa>
> tel : +221 33 827 55 16
> portable: +221 77 123 81 69
>
>
>
>
--
Ahmadou H. DICKO
statistician economist (Ingénieur Statisticien Économiste)
PhD candidate in Climate change economics
Faculty of economics and managment - Cheikh Anta Diop University
West African Science Service Center on Climate Change and Adaptated Land
Use (WASCAL)
Center for Development Research (ZEF) - University of Bonn
email : ahmadou.dicko at ucad.edu.sn
twitter : @dickoah
github : github/dickoa <https://github.com/dickoa>
tel : +221 33 827 55 16
portable: +221 77 123 81 69
[[alternative HTML version deleted]]
More information about the R-sig-Geo
mailing list