[gdal-dev] assigning vertical datum to image files

Even Rouault even.rouault at spatialys.com
Tue Jun 9 07:51:58 PDT 2015


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


More information about the gdal-dev mailing list