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
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
