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

Vinit Bitla bitlavinit at gmail.com
Thu Aug 28 19:43:46 CEST 2014 

Hello All,

Can any one can help me to convert this Lidar data latitude and longitude
in Degree Decimal format using R.

  latitude longitude elevation intensity

1 204140.7   4733835    334.14       6.4
2 204150.1   4733849    334.86      11.4

Appreciate your help.


Vinit



On Thu, Aug 14, 2014 at 3: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
>
>   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
>

	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list