naturalatlas / node-gdal

Node.js bindings for GDAL (Geospatial Data Abstraction Library)
http://naturalatlas.github.io/node-gdal/classes/gdal.html
Apache License 2.0
566 stars 124 forks source link

Example "gdalinfo" implementation gives different corner coordinates #206

Open perhallstroem opened 7 years ago

perhallstroem commented 7 years ago

Hi, I have noticed that the coordinates that I get with the example implementation of "gdalinfo" is different from those that my installed gdalinfo binary/command gives me.

Relevant sections of the output:

From "real" gdalinfo:

Upper Left  (  446246.164,  336322.915) (  1d18'43.42"W, 52d55'19.24"N)
Lower Left  (  446246.164,  329294.211) (  1d18'47.02"W, 52d51'31.77"N)
Upper Right (  447593.502,  336322.915) (  1d17'31.28"W, 52d55'18.81"N)
Lower Right (  447593.502,  329294.211) (  1d17'34.98"W, 52d51'31.34"N)
Center      (  446919.833,  332808.563) (  1d18' 9.18"W, 52d53'25.29"N)

From example gdalinfo:

Upper Left   (446246.16, 336322.91) (  1d18'48.92"W,  52d55'20.41"N)
Upper Right  (447593.5, 336322.91) (  1d17'36.78"W,  52d55'19.99"N)
Bottom Right (447593.5, 329294.21) (  1d17'40.48"W,  52d51'32.55"N)
Bottom Left  (446246.16, 329294.21) (  1d18'52.51"W,  52d51'32.97"N)
Center       (446919.83, 332808.56) (  1d18'14.67"W,  52d53'26.48"N)

Notice that the coordinates in the native CRS (EPSG:27700) are the same, but the transformed WGS84 coordinates differ by quite a substantial amount. Not sure what other information you need. Please let me know and I'll attempt to be as helpful as possible.

WKT:

PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646000044,
                AUTHORITY["EPSG","7001"]],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27700"]]
perhallstroem commented 7 years ago

Okay, so after reading a bit more code it is clear that gdalinfo actually does not print WGS84 coordinates but coordinates based on the GEOGCS of the projected coordinate system. To get the same numbers, this function roughly mimics the behaviour of gdalinfo, as far as I can understand:

const getBestGeoCs = (ds) => {
  const wgs84 = gdal.SpatialReference.fromEPSG(4326);
  if(!ds.srs) return wgs84;

  const cloneGeogCS = ds.srs.cloneGeogCS();
  return !cloneGeogCS ? wgs84 : cloneGeogCS;
};