GDAL's coordinate system demo

osr_test.c

#include <ogr_srs_api.h>
#include <cpl_conv.h>

int main(void) {
    OGRSpatialReferenceH sr = OSRNewSpatialReference(NULL);
    OSRImportFromEPSG(sr, 3857);
//    OSRSetWellKnownGeogCS(sr, "WGS84");
    char* ppsz;
    OSRExportToPrettyWkt(sr, &ppsz, FALSE);
    printf("3857.wkt: \n%s\n", ppsz);
    CPLFree(ppsz);
    printf("3857.au: %s\n", OSRGetAuthorityCode(sr,NULL));
    printf("3857.geogcs.au: %s\n", OSRGetAuthorityCode(sr, "GEOGCS"));
    printf("3857.geogcs|unit.au: %s\n", OSRGetAuthorityCode(sr, "GEOGCS|UNIT"));
    printf("3857.geogcs|datum.au: %s\n", OSRGetAuthorityCode(sr, "GEOGCS|DATUM"));
    printf("3857.geogcs|datum|spheroid.au: %s\n", OSRGetAuthorityCode(sr, "GEOGCS|DATUM|SPHEROID"));
    printf("3857.projcs.au: %s\n", OSRGetAuthorityCode(sr, "PROJCS"));

    //OGRSpatialReferenceH hLatLong = OSRCloneGeogCS(sr);
    OGRSpatialReferenceH hLatLong = OSRNewSpatialReference(NULL);
    OSRSetWellKnownGeogCS(hLatLong, "WGS84");
    printf("wgs84.au: %s\n", OSRGetAuthorityCode(hLatLong, NULL));
    printf("wgs84.geogcs.au: %s\n", OSRGetAuthorityCode(hLatLong, "GEOGCS"));
    printf("wgs84.geogcs|unit.au: %s\n", OSRGetAuthorityCode(hLatLong, "GEOGCS|UNIT"));
    printf("wgs84.geogcs|datum.au: %s\n", OSRGetAuthorityCode(hLatLong, "GEOGCS|DATUM"));
    printf("wgs84.geogcs|datum|spheroid.au: %s\n", OSRGetAuthorityCode(hLatLong, "GEOGCS|DATUM|SPHEROID"));
    printf("wgs84.projcs.au: %s\n", OSRGetAuthorityCode(hLatLong, "PROJCS"));

    OGRCoordinateTransformationH hTransform = OCTNewCoordinateTransformation(sr,hLatLong);

    double x = 12955467.346012;
    double y = 4854848.34787838;
    double z = 0;
    OCTTransform(hTransform, 1, &x, &y, &z);

    printf("x: %f, y: %f, z: %f\n", x, y, z);
    OCTDestroyCoordinateTransformation(hTransform);
    OSRDestroySpatialReference(sr);

    OGRSpatialReferenceH hUTM = OSRNewSpatialReference(NULL);
    OSRSetProjCS(hUTM, "UTM 20 / WGS84");
    OSRSetWellKnownGeogCS(hUTM, "WGS84");
    OSRSetUTM(hUTM, 20, TRUE);

    hTransform = OCTNewCoordinateTransformation(hLatLong, hUTM);
    OCTTransform(hTransform, 1, &x, &y, &z);
    printf("x: %f, y: %f, z: %f\n", x, y, z);

    OCTDestroyCoordinateTransformation(hTransform);
    OSRDestroySpatialReference(hUTM);
    
    return 0;
}

gcc gdal-config --cflags osr_test.c gdal-config --libs -o osr_test