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.

_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

Reply via email to