[gdal-dev] assigning vertical datum to image files
Newcomb, Doug
doug_newcomb at fws.gov
Tue Jun 9 12:07:34 PDT 2015
putting the gtx files in the right proj directory works wonders :-)
Thanks! All is working now.
Doug
On Tue, Jun 9, 2015 at 10:51 AM, Even Rouault <even.rouault at spatialys.com>
wrote:
> Selon "Newcomb, Doug" <doug_newcomb at fws.gov>:
>
> > Evan,
> > Here is the output from gdalinfo wit the expanded config.
> > Note the error messages on failure to load datum shift files for the
> > corner coordinates.
> > I'll upgrade to RC1 and check that again, but thought I would mention.
>
> That's not related to GDAL but to PROJ.4. You'd need to download and
> install
> those grid files in the PROJ_LIB directory.
>
> Anyway, I realize I gave you the code for a vertical CS with unit = meters
> whereas yours was ft. So you'd better use 6360 instead (which has no
> reference
> to grids)
>
> >
> > gdalinfo --config GTIFF_REPORT_COMPD_CS YES D05_37_30081003_20141231.tif
> >
> > Driver: GTiff/GeoTIFF
> > Files: D05_37_30081003_20141231.tif
> > Size is 1000, 1000
> > Coordinate System is:
> > COMPD_CS["unknown",
> > PROJCS["NAD83(2011) / North Carolina (ftUS)",
> > GEOGCS["NAD83(2011)",
> > DATUM["NAD83_National_Spatial_Reference_System_2011",
> > SPHEROID["GRS 1980",6378137,298.257222101,
> > AUTHORITY["EPSG","7019"]],
> > AUTHORITY["EPSG","1116"]],
> > PRIMEM["Greenwich",0,
> > AUTHORITY["EPSG","8901"]],
> > UNIT["degree",0.0174532925199433,
> > AUTHORITY["EPSG","9122"]],
> > AUTHORITY["EPSG","6318"]],
> > PROJECTION["Lambert_Conformal_Conic_2SP"],
> > PARAMETER["standard_parallel_1",36.16666666666666],
> > PARAMETER["standard_parallel_2",34.33333333333334],
> > PARAMETER["latitude_of_origin",33.75],
> > PARAMETER["central_meridian",-79],
> > PARAMETER["false_easting",2000000],
> > PARAMETER["false_northing",0],
> > UNIT["US survey foot",0.3048006096012192,
> > AUTHORITY["EPSG","9003"]],
> > AXIS["X",EAST],
> > AXIS["Y",NORTH],
> > AUTHORITY["EPSG","6543"]],
> > VERT_CS["NAVD88 height",
> > VERT_DATUM["North American Vertical Datum 1988",2005,
> >
> >
>
> EXTENSION["PROJ4_GRIDS","g2012a_conus.gtx,g2012a_alaska.gtx,g2012a_guam.gtx,g2012a_hawaii.gtx,g2012a_puertorico.gtx,g2012a_samoa.gtx"],
> > AUTHORITY["EPSG","5103"]],
> > UNIT["metre",1,
> > AUTHORITY["EPSG","9001"]],
> > AXIS["Up",UP],
> > AUTHORITY["EPSG","5703"]]]
> > Origin = (3010000.000000000000000,805000.000000000000000)
> > Pixel Size = (5.000000000000000,-5.000000000000000)
> > Metadata:
> > AREA_OR_POINT=Area
> > DataType=Generic
> > Image Structure Metadata:
> > COMPRESSION=DEFLATE
> > INTERLEAVE=BAND
> > Corner Coordinates:
> > ERROR 1: failed to load datum shift file
> > Upper Left ( 3010000.000, 805000.000)
> > ERROR 1: failed to load datum shift file
> > Lower Left ( 3010000.000, 800000.000)
> > ERROR 1: failed to load datum shift file
> > Upper Right ( 3015000.000, 805000.000)
> > ERROR 1: failed to load datum shift file
> > Lower Right ( 3015000.000, 800000.000)
> > ERROR 1: failed to load datum shift file
> > Center ( 3012500.000, 802500.000)
> > Band 1 Block=1000x2 Type=Float32, ColorInterp=Gray
> > Description = Layer_1
> > Min=-3.000 Max=-3.000
> > Minimum=-3.000, Maximum=-3.000, Mean=-3.000, StdDev=0.000
> > NoData Value=-3.40282306073709653e+38
> > Metadata:
> > LAYER_TYPE=athematic
> > STATISTICS_COVARIANCES=0
> > STATISTICS_MAXIMUM=-3
> > STATISTICS_MEAN=-3
> > STATISTICS_MEDIAN=0
> > STATISTICS_MINIMUM=-3
> > STATISTICS_MODE=0
> > STATISTICS_SKIPFACTORX=1
> > STATISTICS_SKIPFACTORY=1
> > STATISTICS_STDDEV=0
> >
> >
> > Doug
> >
> > On Tue, Jun 9, 2015 at 10:21 AM, Newcomb, Doug <doug_newcomb at fws.gov>
> wrote:
> >
> > > Evan,
> > > I went back and checked. The .aux.xml sidecar files with the img
> files
> > > do not appear to have projection information.
> > >
> > > Doug
> > >
> > > On Tue, Jun 9, 2015 at 9:32 AM, Even Rouault <
> even.rouault at spatialys.com>
> > > wrote:
> > >
> > >> Selon "Newcomb, Doug" <doug_newcomb at fws.gov>:
> > >>
> > >> > Thanks!
> > >> >
> > >> > That should work in this case to build the virtual raster. I'm
> still
> > >> > curious about assigning the vertical datum to the rasters.
> > >>
> > >> Doug,
> > >>
> > >> I'm a bit surprised that the HFA driver reports a SRS with a VERTCS
> node
> > >> as a
> > >> child of the PROJCS node. I believe this is an incorrect WKT
> definition.
> > >> There
> > >> should rather be a COMPD_CS root note with 2 children: a PROJCS node
> for
> > >> the
> > >> horizontal CRS and a VERTCS node for the vertical CRS. Perhaps this
> SRS
> > >> string
> > >> is defined in the .aux.xml sidecar files that are apparently present ?
> > >>
> > >> Regarding your question, there are a few EPSG codes that correspond to
> > >> aCOMPD_CS, for example 6349. Otherwise, GDAL allows to build a custom
> > >> COMPD_CS
> > >> with the syntax EPSG:XXXX+YYYY where XXXX is the EPSG code of the
> > >> horizontal SRS
> > >> and YYYY the EPSG code of the vertical SRS. In your case that should
> be
> > >> EPSG:6543+5703 I think
> > >>
> > >> So you can do things like "gdal_translate in.tif out.tif -a_srs
> > >> EPSG:6543+5703"
> > >>
> > >> If you read back out.tif, you will only get the horizontal SRS by
> > >> default. This
> > >> was for backward compatibility reason I think. But if you specify
> > >> "--config
> > >> GTIFF_REPORT_COMPD_CS YES" you will get the COMPD_CS.
> > >>
> > >> Even
> > >>
> > >> >
> > >> > Doug
> > >> >
> > >> > On Tue, Jun 9, 2015 at 1:48 AM, Hermann Peifer <peifer at gmx.eu>
> wrote:
> > >> >
> > >> > >
> > >> > > About your initial problem:
> > >> > > > but gdalbuildvrt complained, informing me that
> > >> > > > they were not in the same projection.
> > >> > >
> > >> > > What about using the below gdalbuildvrt option?
> > >> > >
> > >> > > Hermann
> > >> > >
> > >> > > -allow_projection_difference:
> > >> > > (starting with GDAL 1.7.0) When this option is specified, the
> utility
> > >> will
> > >> > > accept to make a VRT even if the input datasets have not the same
> > >> > > projection. Note: this does not mean that they will be
> reprojected.
> > >> Their
> > >> > > projection will just be ignored.
> > >> > >
> > >> > > Source: http://www.gdal.org/gdalbuildvrt.html
> > >> > >
> > >> > >
> > >> > >
> > >> > > On 2015-06-08 21:01, Newcomb, Doug wrote:
> > >> > >
> > >> > >> Hi Folks,
> > >> > >> I have a directory of 783 .img format Lidar - based DEMs in the
> same
> > >> > >> projection, North Carolina State Plane Feet (2011) , NAD 83 ,
> > >> NVD88. I
> > >> > >> was going to use gdalbuildvrt to create a single virtual image
> for
> > >> the
> > >> > >> area, but gdalbuildvrt complained, informing me that they were
> not in
> > >> > >> the same projection.
> > >> > >>
> > >> > >> A couple of quick bash scripts/commands
> > >> > >>
> > >> > >> for x in *.img; do gdalinfo $x >>img_info.txt;done
> > >> > >>
> > >> > >> and
> > >> > >> grep PROJCS img_info.txt|sort|uniq -c
> > >> > >>
> > >> > >> gives me the following:
> > >> > >>
> > >> > >> 437
> > >> PROJCS["NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US",
> > >> > >> 346
> > >> PROJCS["NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet_2011",
> > >> > >>
> > >> > >> gdalinfo gives the following for each type of file ( origin
> section
> > >> > >> clipped out) :
> > >> > >>
> > >> > >> Driver: HFA/Erdas Imagine Images (.img)
> > >> > >> Files: D05_37_20878102_20141231.img
> > >> > >> D05_37_20878102_20141231.img.aux.xml
> > >> > >> Size is 1000, 1000
> > >> > >> Coordinate System is:
> > >> > >> PROJCS["NAD_1983_2011_StatePlane_North_Carolina_FIPS_3200_Ft_US",
> > >> > >> GEOGCS["GCS_NAD_1983_2011",
> > >> > >> DATUM["NAD_1983_2011",
> > >> > >> SPHEROID["GRS_1980",6378137.0,298.257222101]],
> > >> > >> PRIMEM["Greenwich",0.0],
> > >> > >> UNIT["Degree",0.017453292519943295]],
> > >> > >> PROJECTION["Lambert_Conformal_Conic_2SP"],
> > >> > >> PARAMETER["False_Easting",2000000.0],
> > >> > >> PARAMETER["False_Northing",0.0],
> > >> > >> PARAMETER["Central_Meridian",-79.0],
> > >> > >> PARAMETER["Standard_Parallel_1",34.3333333333333],
> > >> > >> PARAMETER["Standard_Parallel_2",36.1666666666667],
> > >> > >> PARAMETER["Latitude_Of_Origin",33.75],
> > >> > >> UNIT["Foot_US",0.30480060960121924],
> > >> > >> VERTCS["NAVD_1988_Foot_US",
> > >> > >> VDATUM["North_American_Vertical_Datum_1988"],
> > >> > >> PARAMETER["Vertical_Shift",0.0],
> > >> > >> PARAMETER["Direction",1.0],
> > >> > >> UNIT["Foot_US",0.3048006096012192]]]
> > >> > >>
> > >> > >>
> > >> > >> Driver: HFA/Erdas Imagine Images (.img)
> > >> > >> Files: D05_37_20957301_20141231.img
> > >> > >> D05_37_20957301_20141231.img.aux.xml
> > >> > >> Size is 1000, 1000
> > >> > >> Coordinate System is:
> > >> > >> PROJCS["NAD_1983_StatePlane_North_Carolina_FIPS_3200_Feet_2011",
> > >> > >> GEOGCS["GCS_NAD_1983_2011",
> > >> > >> DATUM["NAD_1983_2011",
> > >> > >> SPHEROID["GRS_1980",6378137.0,298.257222101]],
> > >> > >> PRIMEM["Greenwich",0.0],
> > >> > >> UNIT["Degree",0.0174532925199433]],
> > >> > >> PROJECTION["Lambert_Conformal_Conic_2SP"],
> > >> > >> PARAMETER["False_Easting",2000000.002616666],
> > >> > >> PARAMETER["False_Northing",0.0],
> > >> > >> PARAMETER["Central_Meridian",-79.0],
> > >> > >> PARAMETER["Standard_Parallel_1",34.33333333333334],
> > >> > >> PARAMETER["Standard_Parallel_2",36.16666666666666],
> > >> > >> PARAMETER["Latitude_Of_Origin",33.75],
> > >> > >> UNIT["Foot_US",0.30480060960121924],
> > >> > >> VERTCS["NAVD_1988_Foot_US",
> > >> > >> VDATUM["North_American_Vertical_Datum_1988"],
> > >> > >> PARAMETER["Vertical_Shift",0.0],
> > >> > >> PARAMETER["Direction",1.0],
> > >> > >> UNIT["Foot_US",0.3048006096012192]]]
> > >> > >>
> > >> > >>
> > >> > >>
> > >> > >> In theory, these should all be EPSG:6543, so just assigning the
> > >> correct
> > >> > >> horizontal projection/datum with the epsg code should make
> things
> > >> > >> usable. (i.e,
> > >> > >> gdal_translate -a_"srs epsg:6543" --config GDAL_CACHEMAX 512 -of
> > >> GTiff
> > >> > >> -co COMPRESS=DEFLATE -co PREDICTOR=3 in.img out.tif ) ( Thank
> you
> > >> for
> > >> > >> the EPSG:6543 projection support in GDAL 2.0!)
> > >> > >>
> > >> > >> However, this only assigns the horizontal infomation, how does
> one
> > >> > >> assign a vertical datum with a horizontal EPSG code?
> > >> > >>
> > >> > >> I see the NVD88 code in vertcs.csv , but how would I implement
> it in
> > >> the
> > >> > >> above command?
> > >> > >>
> > >> > >> Using gdal 2.0 Beta2
> > >> > >>
> > >> > >>
> > >> > >> Doug
> > >> > >> --
> > >> > >> Doug Newcomb
> > >> > >> USFWS
> > >> > >> Raleigh, NC
> > >> > >> 919-856-4520 ext. 14 doug_newcomb at fws.gov <mailto:
> > >> doug_newcomb at fws.gov>
> > >> > >>
> > >> > >>
> > >> >
> > >>
> > >>
> >
>
> ---------------------------------------------------------------------------------------------------------
> > >> > >> The opinions I express are my own and are not representative of
> the
> > >> > >> official policy of the U.S.Fish and Wildlife Service or Dept. of
> the
> > >> > >> Interior. Life is too short for undocumented, proprietary data
> > >> formats.
> > >> > >>
> > >> > >>
> > >> > >> _______________________________________________
> > >> > >> gdal-dev mailing list
> > >> > >> gdal-dev at lists.osgeo.org
> > >> > >> http://lists.osgeo.org/mailman/listinfo/gdal-dev
> > >> > >>
> > >> > >>
> > >> > >
> > >> >
> > >> >
> > >> > --
> > >> > Doug Newcomb
> > >> > USFWS
> > >> > Raleigh, NC
> > >> > 919-856-4520 ext. 14 doug_newcomb at fws.gov
> > >> >
> > >>
> > >>
> >
>
> ---------------------------------------------------------------------------------------------------------
> > >> > The opinions I express are my own and are not representative of the
> > >> > official policy of the U.S.Fish and Wildlife Service or Dept. of the
> > >> > Interior. Life is too short for undocumented, proprietary data
> > >> formats.
> > >> >
> > >>
> > >>
> > >> --
> > >> Spatialys - Geospatial professional services
> > >> http://www.spatialys.com
> > >>
> > >
> > >
> > >
> > > --
> > > Doug Newcomb
> > > USFWS
> > > Raleigh, NC
> > > 919-856-4520 ext. 14 doug_newcomb at fws.gov
> > >
> > >
> >
>
> ---------------------------------------------------------------------------------------------------------
> > > The opinions I express are my own and are not representative of the
> > > official policy of the U.S.Fish and Wildlife Service or Dept. of the
> > > Interior. Life is too short for undocumented, proprietary data
> formats.
> > >
> >
> >
> >
> > --
> > Doug Newcomb
> > USFWS
> > Raleigh, NC
> > 919-856-4520 ext. 14 doug_newcomb at fws.gov
> >
>
> ---------------------------------------------------------------------------------------------------------
> > The opinions I express are my own and are not representative of the
> > official policy of the U.S.Fish and Wildlife Service or Dept. of the
> > Interior. Life is too short for undocumented, proprietary data formats.
> >
>
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
>
--
Doug Newcomb
USFWS
Raleigh, NC
919-856-4520 ext. 14 doug_newcomb at fws.gov
---------------------------------------------------------------------------------------------------------
The opinions I express are my own and are not representative of the
official policy of the U.S.Fish and Wildlife Service or Dept. of the
Interior. Life is too short for undocumented, proprietary data formats.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20150609/b7c4578b/attachment-0001.html>
More information about the gdal-dev
mailing list