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

[Edzer Pebesma]

Thanks, Even!

I was after the units of the horizontal coordinates; GetLinearUnits() 
returns "unknown" for EPSG:4326 and, according to the docs returns the 
units of the vertical coordinate in 3D CRS (didn't check). So what I now 
use is GetAttrValue("UNIT", 0) (which returns "degree" for EPSG:4326"), 
and if that returns NULL, use GetLinearUnit(&unit) if its value differs 
from "unknown".

Many regards, and thanks for the example code,

On 07/03/2023 12:46, Even Rouault wrote:
> 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?
> 

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081