mmomtchev / node-gdal-async

Node.js bindings for GDAL (Geospatial Data Abstraction Library) with full async support
https://mmomtchev.github.io/node-gdal-async/
Apache License 2.0
129 stars 26 forks source link

extend serve.js to show WGS84 transform? #70

Closed Pomax closed 1 year ago

Pomax commented 1 year ago

https://github.com/mmomtchev/node-gdal-async/blob/master/examples/serve.js currently shows how to get pixels from a geotiff, but does not show how to do so while taking a standard projection like WGS84 into account. Could the example be updated to show how an incoming lat/long gets mapped to its corresponding x/y pixel array position with a coordinate transform taken into account?

Pomax commented 1 year ago

awesome, thanks @mmomtchev!

Pomax commented 1 year ago

Hmm, something might not be quite right: if you use that with an ALOS World3D tif (specifically, ALPSMLC30_S045E168_DSM.tif) and I curl http://localhost:8010/WGS84/-44.9678876/168.3115521, I get

Error: ALPSMLC30_S045E168_DSM.tif, band 1: Access window out of range in RasterIO().  Requested
(-766685,-764322) of size 1x1 on raster of 3600x3600.
Pomax commented 1 year ago

A quick info check confirms this TIF uses WGS84:

>node gdaltest.js
Opened cache\ALPSMLC30_S045E168_DSM.tif
Coordinate System is: 
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]
Origin = (168, -44)
Pixel Size = (0.0002777777777777778, -0.0002777777777777778)
GeoTransform =
[ 168, 0.0002777777777777778, 0, -44, 0, -0.0002777777777777778 ]
Pomax commented 1 year ago

Ahhhhh lat/long reversal! It looks like where the example currently says const coords = xform.transformPoint(+params[2], +params[3]); it should say const coords = xform.transformPoint(+params[3], +params[2]); instead. Is there a way to check which ordering a file needs?

mmomtchev commented 1 year ago

With EPSG:4326, transformPoint always expect x, y order which it transforms to the correct x, y coordinates of the Dataset.