Closed graemejcox closed 4 years ago
Hi @graemejcox!
Thanks for your interest in geotiff.js.
I'm not able to provide a full answer to this question, as "CanvasLayer.Field" is not part of this project, so you have to find help elsewhere (I think you are refering to this project: https://github.com/IHCantabria/Leaflet.CanvasLayer.Field).
On the side of geotiff.js I can say that the full pixel transform matrix (i.e GeoTransform
in GDAL) can be extracted and inspected. Everything else depends on the rendering engine.
Sorry, I meant Geotiff.js not CanvasLayer.Field. Good to see Geotiff.js supports transform matrix. Thanks of your efforts.
I have some geotiffs that have rotation or skew. Can CanvasLayer.Field handle these?
This is how we set it with GDAL SatGeoTransform: Given information from the aforementioned gdal datamodel docs, the 3rd & 5th parameters of SatGeoTransform (x_skew and y_skew respectively) can be calculated from two control points (p1, p2) with known x and y in both "geo" and "pixel" coordinate spaces. p1 should be above-left of p2 in pixelspace.
x_skew = sqrt((p1.geox-p2.geox)2 + (p1.geoy-p2.geoy)2) / (p1.pixely - p2.pixely) y_skew = sqrt((p1.geox-p2.geox)2 + (p1.geoy-p2.geoy)2) / (p1.pixelx - p2.pixelx) In short this is the ratio of Euclidean distance between the points in geospace to the height (or width) of the image in pixelspace.
The units of the parameters are "geo"length/"pixel"length.
Here is a demonstration using the corners of the image stored as control points (gcps):
import gdal from math import sqrt
ds = gdal.Open(fpath) gcps = ds.GetGCPs()
assert gcps[0].Id == 'UpperLeft' p1 = gcps[0] assert gcps[2].Id == 'LowerRight' p2 = gcps[2]
y_skew = ( sqrt((p1.GCPX-p2.GCPX)2 + (p1.GCPY-p2.GCPY)2) / (p1.GCPPixel - p2.GCPPixel) ) x_skew = ( sqrt((p1.GCPX-p2.GCPX)2 + (p1.GCPY-p2.GCPY)2) / (p1.GCPLine - p2.GCPLine) )
x_res = (p2.GCPX - p1.GCPX) / ds.RasterXSize y_res = (p2.GCPY - p1.GCPY) / ds.RasterYSize
ds.SetGeoTransform([ p1.GCPX, x_res, x_skew, p1.GCPY, y_skew, y_res, ])