[gdal-dev] getting SRS units from GetAttrValue("UNIT", 0)?

[Even Rouault]

Edzer,

if you change your code to use GetLinearUnits() it will work:

#include <gdal.h>
#include <gdal_priv.h>
#include <iostream>

int main()
{
     OGRSpatialReference *srs = new OGRSpatialReference;
     int epsg = 3031;
     srs->importFromEPSG(epsg);
     const char* unit = "";
     srs->GetLinearUnits(&unit);
     std::cout << epsg << ":" << unit << std::endl;
     const char *proj4s = "+proj=eqearth";
     srs->importFromProj4(proj4s);
     unit = "";
     srs->GetLinearUnits(&unit);
     std::cout << proj4s << ":" << unit << std::endl;
     return 0;
}


If you wonder why GetAttrValue("UNIT", 0) works sometimes and not not 
for others, this is a bit tricky. When querying WKT nodes directly like 
with GetAttrValue(), GDAL works internally with a WKT1 representation, 
for backward compatibility as a lot of code inside GDAL (and potentially 
outside) queries with WKT1 node names. The issue here is that eqearth is 
a new projection method that didn't exist in GDAL <= 2.4, so GDAL 
fallbacks to WKT2 representation, and it's not the UNIT keyword in WKT2 
but LENGTHUNIT.

Using GetLinearUnits() is more robust, as independent from the WKT1/WKT2 
representation

Even


Le 07/03/2023 à 11:38, Edzer Pebesma a écrit :
> I switched getting coordinate units from the expanded proj4string 
> representation to GetAttrValue("UNIT", 0), but failed; try
>
> #include <gdal.h>
> #include <gdal_priv.h>
> #include <iostream>
>
> int main()
> {
>     OGRSpatialReference *srs = new OGRSpatialReference;
>     int epsg = 3031;
>     srs->importFromEPSG(epsg);
>     std::cout << epsg << ":" << srs->GetAttrValue("UNIT",0) << std::endl;
>     const char *proj4s = "+proj=eqearth";
>     srs->importFromProj4(proj4s);
>     const char *u = srs->GetAttrValue("UNIT",0);
>     std::cout << proj4s << ":" << (u ? u : "NULL") << std::endl;
>     return 0;
> }
>
> which gives
>
> 3031:metre
> +proj=eqearth:NULL
>
> whereas the proj4string expansion of "+proj=eqearth" is "+proj=eqearth 
> +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs", so has the units.
>
> What do I miss here? What is the recommended way to get coordinate units?

-- 
http://www.spatialys.com
My software is free, but my time generally not.