5

In Bash scripts, I use gdal(srs)info / ogrinfo to get information from raster / vector files.

To use this information in my scripts, I'm extracting information using grep

... e.g. like this (to receive ulx)

gdalinfo raster.tif | grep "Upper Left" | grep -Po '[-0-9\.]+(?=,)'

... or this (to receive EPSG Code)

gdalsrsinfo raster.tif | grep -Po '(?<=AUTHORITY\["EPSG",")\d+(?="\]\]$)'

That works, but I'm not very happy with that solution because it's very dependent on the way gdalinfo works. If anything changes, my scripts might break.

Does anyone know of any better solution to that ?

whyzar
  • 12,053
  • 23
  • 38
  • 72
pLumo
  • 6,489
  • 20
  • 41
  • Change from bash to a script language with GDAL/ORG bindings (for example Python). Here is an example: http://gis.stackexchange.com/questions/57834/how-to-get-raster-corner-coordinates-using-python-gdal-bindings – Zoltan Jan 03 '17 at 15:18
  • Perhaps rio info from rasterio cli can help. rio info raster.tif --crs will directly give you the proj4 string for instance. See rio info --help – Loïc Dutrieux Jan 04 '17 at 00:49

3 Answers3

5

gdalsrsinfo has an output format option (-o) which has a range of formats. epsg is still experimental as of Jan 17, and I'd not recommend it at this stage. xml will probably tell you more than you possibly wanted to know:

$ gdalsrsinfo -o xml phoenix_DEM_subset.tif 
<gml:GeographicCRS gml:id="ogrcrs1">
  <gml:srsName>WGS 84</gml:srsName>
  <gml:srsID>
    <gml:name codeSpace="urn:ogc:def:crs:EPSG::">4326</gml:name>
  </gml:srsID>
  <gml:usesEllipsoidalCS>
    <gml:EllipsoidalCS gml:id="ogrcrs2">
      <gml:csName>ellipsoidal</gml:csName>
      <gml:csID>
        <gml:name codeSpace="urn:ogc:def:cs:EPSG::">6402</gml:name>
      </gml:csID>
      <gml:usesAxis>
        <gml:CoordinateSystemAxis gml:id="ogrcrs3" gml:uom="urn:ogc:def:uom:EPSG::9102">
          <gml:name>Geodetic latitude</gml:name>
          <gml:axisID>
            <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9901</gml:name>
          </gml:axisID>
          <gml:axisAbbrev>Lat</gml:axisAbbrev>
          <gml:axisDirection>north</gml:axisDirection>
        </gml:CoordinateSystemAxis>
      </gml:usesAxis>
      <gml:usesAxis>
        <gml:CoordinateSystemAxis gml:id="ogrcrs4" gml:uom="urn:ogc:def:uom:EPSG::9102">
          <gml:name>Geodetic longitude</gml:name>
          <gml:axisID>
            <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9902</gml:name>
          </gml:axisID>
          <gml:axisAbbrev>Lon</gml:axisAbbrev>
          <gml:axisDirection>east</gml:axisDirection>
        </gml:CoordinateSystemAxis>
      </gml:usesAxis>
    </gml:EllipsoidalCS>
  </gml:usesEllipsoidalCS>
  <gml:usesGeodeticDatum>
    <gml:GeodeticDatum gml:id="ogrcrs5">
      <gml:datumName>WGS_1984</gml:datumName>
      <gml:datumID>
        <gml:name codeSpace="urn:ogc:def:datum:EPSG::">6326</gml:name>
      </gml:datumID>
      <gml:usesPrimeMeridian>
        <gml:PrimeMeridian gml:id="ogrcrs6">
          <gml:meridianName>Greenwich</gml:meridianName>
          <gml:greenwichLongitude>
            <gml:angle uom="urn:ogc:def:uom:EPSG::9102">0</gml:angle>
          </gml:greenwichLongitude>
        </gml:PrimeMeridian>
      </gml:usesPrimeMeridian>
      <gml:usesEllipsoid>
        <gml:Ellipsoid gml:id="ogrcrs7">
          <gml:ellipsoidName>WGS 84</gml:ellipsoidName>
          <gml:ellipsoidID>
            <gml:name codeSpace="urn:ogc:def:ellipsoid:EPSG::">7030</gml:name>
          </gml:ellipsoidID>
          <gml:semiMajorAxis uom="urn:ogc:def:uom:EPSG::9001">6378137</gml:semiMajorAxis>
          <gml:secondDefiningParameter>
            <gml:inverseFlattening uom="urn:ogc:def:uom:EPSG::9201">298.257223563</gml:inverseFlattening>
          </gml:secondDefiningParameter>
        </gml:Ellipsoid>
      </gml:usesEllipsoid>
    </gml:GeodeticDatum>
  </gml:usesGeodeticDatum>
</gml:GeographicCRS>

If you have a more recent version of gdalinfo (I tested on master from Dec 16, but anything from 2.0 onwards should be OK), you can use the -json format option to get verbose descriptive output.

$ gdalinfo -json phoenix_DEM_subset.tif 
{
  "description":"phoenix_DEM_subset.tif",
  "driverShortName":"GTiff",
  "driverLongName":"GeoTIFF",
  "files":[
    "phoenix_DEM_subset.tif"
  ],
  "size":[
    4633,
    1153
  ],
  "coordinateSystem":{
    "wkt":"GEOGCS[\"WGS 84\",\n    DATUM[\"WGS_1984\",\n        SPHEROID[\"WGS 84\",6378137,298.257223563,\n            AUTHORITY[\"EPSG\",\"7030\"]],\n        AUTHORITY[\"EPSG\",\"6326\"]],\n    PRIMEM[\"Greenwich\",0],\n    UNIT[\"degree\",0.0174532925199433],\n    AUTHORITY[\"EPSG\",\"4326\"]]"
  },
  "geoTransform":[
    -112.114887330000002,
    0.000083333333333,
    0.0,
    33.4441373299999967,
    0.0,
    -0.000083333333333
  ],
  "metadata":{
    "":{
      "AREA_OR_POINT":"Area",
      "TIFFTAG_XRESOLUTION":"1",
      "TIFFTAG_YRESOLUTION":"1"
    },
    "IMAGE_STRUCTURE":{
      "INTERLEAVE":"BAND"
    }
  },
  "cornerCoordinates":{
    "upperLeft":[
      -112.1148873,
      33.4441373
    ],
    "lowerLeft":[
      -112.1148873,
      33.348054
    ],
    "upperRight":[
      -111.728804,
      33.4441373
    ],
    "lowerRight":[
      -111.728804,
      33.348054
    ],
    "center":[
      -111.9218457,
      33.3960957
    ]
  },
  "wgs84Extent":{
    "type":"Polygon",
    "coordinates":[
      [
        [
          -112.1148873,
          33.4441373
        ],
        [
          -112.1148873,
          33.348054
        ],
        [
          -111.728804,
          33.4441373
        ],
        [
          -111.728804,
          33.348054
        ],
        [
          -112.1148873,
          33.4441373
        ]
      ]
    ]
  },
  "bands":[
    {
      "band":1,
      "block":[
        4633,
        1
      ],
      "type":"Int16",
      "colorInterpretation":"Gray",
      "metadata":{
      }
    }
  ]

Note that this is not in the 1.x series.

BradHards
  • 12,881
  • 2
  • 37
  • 70
  • I like the json output. =) But then I need to extract the information from the json. That could be done in jq, but I don't have that installed in the production system. Python is installed though. I can use it like this to get ulx: gdalinfo raster.tif -json | python -c "import json,sys; jdata = sys.stdin.read(); data = json.loads(jdata); print(data['cornerCoordinates']['upperLeft'][0]);" But then I could use Python anyways ... – pLumo Jan 04 '17 at 08:13
  • 1
    I'm not sure I see another option here. Either you parse the structured output in a language of your choice, or you use the text format and assume it doesn't change format (which I think it is unlikely to). – BradHards Jan 04 '17 at 08:37
  • You're right. You should add the options to the answer for that I can accept it ;-) But your assumption, that the (json) format is unlikely to change is probably not correct --> for example, but that will affect any means of interpreting – pLumo Jan 04 '17 at 12:36
4

Newer versions of GDAL now have a -json output which can easily be parsed with jq. For example

gdalinfo -json file.tif | jq -r .cornerCoordinates.upperLeft
glaroc
  • 403
  • 2
  • 8
0

This is something I wondered about too.

gdalinfo spits out text to stdout, but in a text format which (as you say) is difficult to parse.

You could use (or port) an existing tool such as gdalinfo-json to convert gdalinfo output into JSON format, and pipe it into a tool such as jq to filter the output.

I haven't tried these, but it might be worth looking at if you don't fancy starting from scratch with python/gdal.

Steven Kay
  • 20,405
  • 5
  • 31
  • 82