Tell me more ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 1 down vote favorite

From GDAL 1.9.0 a new utility is available: gdalsrsinfo.

Although it can yield an XML (GML-based) definition from a source equivalent format or file, it seems that it can only read usual GDAL/OGR forms (complete WKT, PROJ.4, EPSG:n or a file containing the SRS), so no GML definitions. Indeed:

$ gdalsrsinfo -o xml '+proj=longlat +a=3396190 +b=3396190 +no_defs'
<gml:GeographicCRS...

but:

$ gdalsrsinfo -o wkt '<gml:GeographicCRS...</gml:GeographicCRS>'
ERROR 1: ERROR - failed to load SRS definition from <gml:GeographicCRS...

Any of you is aware of a way I could make this translation?

This question is somehow linked to: Can GeoTools parse GML CRS definitions?

share|improve this question

3 Answers

up vote 2 down vote
+50

You are close. Using a modified version of your example works just fine (However, see my updates below to understand why this is happening):

gdalsrsinfo -o xml '+proj=longlat +a=3396190 +b=3396190 +no_defs +rf=1' > /tmp/test1;gdalsrsinfo -o wkt /tmp/test1

I'm taking the output (in GML) of the first command and sending it to a temporary file, then reading it back in wkt format. This demonstrates that what you are asking is indeed possible with gdalsrsinfo.

There is an issue with directly attempting to pipe the output from one gdalsrsinfo call to another because of the way it detects SRS information from a string or a file. May be a bug, I don't know.

Also notice I included +datum=WGS84 to give it something to find in the srs tags. This may be causing you an issue if you were using gml generated from your original test string.

UPDATE: Your issue is being caused because of an inverse flattening value of 0. Add +rf=1 and tweak to your preferred ellipsoid and you will find that this fixes the error. (I'm using 1 as a placeholder of course). I've updated the command string above.

To get an appropriate value, use the formula

f=(a-b)/a which in this case is 0 for you, but this doesn't compute for gdalsrsinfo.

UPDATE 2: If you change your ellipsoid slightly, this will fix the issue as well without having to fake the inverse flattening value. Here's an updated command string that does what you want:

gdalsrsinfo -o xml '+proj=longlat +a=3396195 +b=3396190 +no_defs' > /tmp/test1;gdalsrsinfo -o wkt /tmp/test1

So we see that it doesn't like perfect spheres!

share|improve this answer
Thanks BuckFilledPlatypus, the problem is that if you define a WGS84 datum then your custom datum (the sphere of radius 3396190m in this case) is ignored and replaced with the WGS84 ellipsoid! I guess I am missing some required parameter in the PROJ4 string that inhibits then the back-translation of the XML definition. I will investigate. Any idea or hint meanwhile? – Campa Jul 25 '12 at 10:15
I've updated my original answer above. Your problem is Inverse Flattening. It is not defined in your original string. Add a '+rf=1' and tweak to your preferred ellipsoid, but you will find that this is causing the error. – BuckFilledPlatypus Jul 25 '12 at 16:04
Uh nice.. but this way the WKT would have a spheroid with inverse flattening of 1, which is not my case. In case of a sphere I would have an infinite inverse flattening, but seemingly in the WKT this would appear as SPHEROID["unnamed",3396190,0]. Looking at the example (which really fits my case) in here, I still have the error... Might be some GDAL bug as well? – Campa Jul 26 '12 at 13:46
Please see my updates (UPDATE 2) on the original answer where I explained in more detail what is going on. 1 is a placeholder to demonstrate what is causing your issue not what value you should place there. I'll paste it here as well. gdalsrsinfo -o xml '+proj=longlat +a=3396195 +b=3396190 +no_defs' > /tmp/test1;gdalsrsinfo -o wkt /tmp/test1 . Here we don't need to use the +rf variable because the first operation generates a correct inverse flattening value if the semimajor and semiminor are not the same. – BuckFilledPlatypus Jul 26 '12 at 14:56
Alternatively you can you a sphere with '+ellips=sphere', however the semimajor and semiminor radii are locked to 'a=6370997.0, b=6370997.0' according to the source code. – BuckFilledPlatypus Jul 26 '12 at 15:05
show 2 more comments
up vote 1 down vote

You should be able to find your SRS on http://spatialreference.org/ and they have it in many different formats. E.g. EPSG:4326

Note: Looks like they're having some CSS issues today, but all features seem to be working.

share|improve this answer
Thanks Vadim. Yes I know about spatialreference.org: that's why I thought there is a way out there to translate GML to some common projection string. I would need to do that for arbitrary definitions and on the fly of a web request. – Campa Jul 23 '12 at 10:02
up vote 1 down vote

One undocumented feature is the use of both tools:

kyle@kyle-workstation:Desktop$ gdalsrsinfo http://spatialreference.org/ref/epsg/4326/

PROJ.4 : '+proj=longlat +datum=WGS84 +no_defs '

OGC WKT :
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.