[R-sig-Geo] Converting .hdf file to raster layer
Roger Bivand
Roger.Bivand at nhh.no
Thu Aug 28 21:08:13 CEST 2014
On Thu, 28 Aug 2014, Koen Hufkens wrote:
> Hi Vinti,
>
> Have a look here:
>
> http://gis.stackexchange.com/questions/20018/how-can-i-convert-data-in-the-form-of-lat-lon-value-into-a-raster-file-using-r
No, the question was how to go from a 3D projected data set to
geographical coordinates when the projection is unknown - this is the
opposite direction to the SO thread. I'm also interpreting "in" to mean
"into", as the values are most likely in metres.
This question has no answer unless the exact projection and datum are
provided. Indeed, saying that the desired operation is to go to
geographical coordinates without specifying the target datum is also
fraught with uncertainties (someone has to assume a datum, probably
messing up the Z value totally).
So the correct first step is to find out the coordinate reference system
specification for the input data. Once that is known, everything else
follows.
See for example:
http://www.asprs.org/a/resources/grids/03-2014-India.pdf
(guessed as the OP did not provide an affiliation, but searching shows a
possible location in Pune). As Cliff Mugnier writes, because of secrecy,
the CRS may be hard to decipher. If you reversed lat/long, maybe this is
UTM, but we don't know the datum. Your data source should know.
Roger
PS. The OP has also hijacked an existing thread for this different
question - please read list instructions and the posting guide.
>
> This should get you on your way.
>
> Cheers,
> K
>
>
> On Thu, Aug 28, 2014 at 1:43 PM, Vinit Bitla <bitlavinit at gmail.com> wrote:
>
>> 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]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list