import Log from "../Log"
/**
* The resulting geotransform is the equivalent to padfGT1 and then padfGT2 being applied to a point.
* @param {Float64Array} padfGT1 the first geotransform, six values.
* @param {Float64Array} padfGT2 the second geotransform, six values.
* @param {Float64Array} padfGTOut the output geotransform, six values, may safely be the same array as padfGT1 or padfGT2.
*/
export function GDALComposeGeoTransforms(padfGT1, padfGT2, padfGTOut) {
// We need to think of the geotransform in a more normal form to do
// the matrix multiple:
//
// __ __
// | gt[1] gt[2] gt[0] |
// | gt[4] gt[5] gt[3] |
// | 0.0 0.0 1.0 |
// -- --
//
// Then we can use normal matrix multiplication to produce the
// composed transformation. I don't actually reform the matrix
// explicitly which is why the following may seem kind of spagettish.
padfGTOut[1] = padfGT2[1] * padfGT1[1] + padfGT2[2] * padfGT1[4]
padfGTOut[2] = padfGT2[1] * padfGT1[2] + padfGT2[2] * padfGT1[5]
padfGTOut[0] =
padfGT2[1] * padfGT1[0] + padfGT2[2] * padfGT1[3] + padfGT2[0] * 1.0
padfGTOut[4] = padfGT2[4] * padfGT1[1] + padfGT2[5] * padfGT1[4]
padfGTOut[5] = padfGT2[4] * padfGT1[2] + padfGT2[5] * padfGT1[5]
padfGTOut[3] =
padfGT2[4] * padfGT1[0] + padfGT2[5] * padfGT1[3] + padfGT2[3] * 1.0
}
/**
* Invert Geotransform.
* @param {Float64Array} gt_in Input geotransform (six doubles - unaltered).
* @param {Float64Array} gt_out Output geotransform (six doubles - updated).
* @returns {Number} TRUE on success or FALSE if the equation is uninvertable.
*/
export function GDALInvGeoTransform(gt_in, gt_out) {
// Special case - no rotation - to avoid computing determinate
// and potential precision issues.
if (
gt_in[2] === 0.0 &&
gt_in[4] === 0.0 &&
gt_in[1] !== 0.0 &&
gt_in[5] !== 0.0
) {
/*X = gt_in[0] + x * gt_in[1]
Y = gt_in[3] + y * gt_in[5]
-->
x = -gt_in[0] / gt_in[1] + (1 / gt_in[1]) * X
y = -gt_in[3] / gt_in[5] + (1 / gt_in[5]) * Y
*/
gt_out[0] = -gt_in[0] / gt_in[1]
gt_out[1] = 1.0 / gt_in[1]
gt_out[2] = 0.0
gt_out[3] = -gt_in[3] / gt_in[5]
gt_out[4] = 0.0
gt_out[5] = 1.0 / gt_in[5]
return 1
}
// Assume a 3rd row that is [1 0 0].
// Compute determinate.
const det = gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4]
const magnitude = Math.max(
Math.max(Math.abs(gt_in[1]), Math.abs(gt_in[2])),
Math.max(Math.abs(gt_in[4]), Math.abs(gt_in[5]))
)
if (Math.abs(det) <= 1e-10 * magnitude * magnitude) return 0
const inv_det = 1.0 / det
// Compute adjoint, and divide by determinate.
gt_out[1] = gt_in[5] * inv_det
gt_out[4] = -gt_in[4] * inv_det
gt_out[2] = -gt_in[2] * inv_det
gt_out[5] = gt_in[1] * inv_det
gt_out[0] = (gt_in[2] * gt_in[3] - gt_in[0] * gt_in[5]) * inv_det
gt_out[3] = (-gt_in[1] * gt_in[3] + gt_in[0] * gt_in[4]) * inv_det
return 1
}
/**
* Applies the following computation, converting a (pixel, line) coordinate into a georeferenced (geo_x, geo_y) location.
* @param {Float64Array} padfGeoTransform Six coefficient GeoTransform to apply.
* @param {Number} dfPixel Input pixel position.
* @param {Number} dfLine Input line position.
* @returns {Object} output location
*/
export function GDALApplyGeoTransform(padfGeoTransform, dfPixel, dfLine) {
const pdfGeoX =
padfGeoTransform[0] +
dfPixel * padfGeoTransform[1] +
dfLine * padfGeoTransform[2]
const pdfGeoY =
padfGeoTransform[3] +
dfPixel * padfGeoTransform[4] +
dfLine * padfGeoTransform[5]
return { pdfGeoX, pdfGeoY }
}
/**
* Ground Control Point
* @typedef GCP
* @type {Object}
* @property {String} pszId Unique identifier, often numeric
* @property {String} pszInfo Informational message or ""
* @property {Number} dfGCPPixel Pixel (x) location of GCP on raster
* @property {Number} dfGCPLine Line (y) location of GCP on raster
* @property {Number} dfGCPX X position of GCP in georeferenced space
* @property {Number} dfGCPY Y position of GCP in georeferenced space
* @property {Number} dfGCPZ Elevation of GCP, or zero if not known
*/
/**
* Given a set of GCPs perform first order fit as a geotransform.
* @param {GCP[]} pasGCPs the list of GCP structures.
* @param {Float64Array} padfGeoTransform the six double array in which the affine geotransformation will be returned.
* @param {Boolean} bApproxOK If FALSE the function will fail if the geotransform is not essentially an exact fit (within 0.25 pixel) for all GCPs.
* @return {Boolean} TRUE on success or FALSE if there aren't enough points to prepare a geotransform, the pointers are ill-determined or if bApproxOK is FALSE and the fit is poor.
*/
export function GDALGCPsToGeoTransform(pasGCPs, padfGeoTransform, bApproxOK) {
let dfPixelThreshold = 0.25
if (!bApproxOK) {
// 这里是在读取配置信息 为方便前端实现逻辑默认开启
bApproxOK = true
dfPixelThreshold = 0.25
// bApproxOK = CPLTestBool(
// CPLGetConfigOption("GDAL_GCPS_TO_GEOTRANSFORM_APPROX_OK", "NO"));
// if (!bApproxOK)
// {
// // coverity[tainted_data]
// dfPixelThreshold = CPLAtof(CPLGetConfigOption(
// "GDAL_GCPS_TO_GEOTRANSFORM_APPROX_THRESHOLD", "0.25"));
// }
}
const nGCPCount = pasGCPs.length
if (nGCPCount < 2) {
return false
}
/* -------------------------------------------------------------------- */
/* Recognise a few special cases. */
/* -------------------------------------------------------------------- */
if (nGCPCount === 2) {
if (
pasGCPs[1].dfGCPPixel == pasGCPs[0].dfGCPPixel ||
pasGCPs[1].dfGCPLine == pasGCPs[0].dfGCPLine
)
return false
padfGeoTransform[1] =
(pasGCPs[1].dfGCPX - pasGCPs[0].dfGCPX) /
(pasGCPs[1].dfGCPPixel - pasGCPs[0].dfGCPPixel)
padfGeoTransform[2] = 0.0
padfGeoTransform[4] = 0.0
padfGeoTransform[5] =
(pasGCPs[1].dfGCPY - pasGCPs[0].dfGCPY) /
(pasGCPs[1].dfGCPLine - pasGCPs[0].dfGCPLine)
padfGeoTransform[0] =
pasGCPs[0].dfGCPX -
pasGCPs[0].dfGCPPixel * padfGeoTransform[1] -
pasGCPs[0].dfGCPLine * padfGeoTransform[2]
padfGeoTransform[3] =
pasGCPs[0].dfGCPY -
pasGCPs[0].dfGCPPixel * padfGeoTransform[4] -
pasGCPs[0].dfGCPLine * padfGeoTransform[5]
return true
}
/* -------------------------------------------------------------------- */
/* Special case of 4 corner coordinates of a non-rotated */
/* image. The points must be in TL-TR-BR-BL order for now. */
/* This case helps avoid some imprecision in the general */
/* calculations. */
/* -------------------------------------------------------------------- */
if (
nGCPCount == 4 &&
pasGCPs[0].dfGCPLine == pasGCPs[1].dfGCPLine &&
pasGCPs[2].dfGCPLine == pasGCPs[3].dfGCPLine &&
pasGCPs[0].dfGCPPixel == pasGCPs[3].dfGCPPixel &&
pasGCPs[1].dfGCPPixel == pasGCPs[2].dfGCPPixel &&
pasGCPs[0].dfGCPLine != pasGCPs[2].dfGCPLine &&
pasGCPs[0].dfGCPPixel != pasGCPs[1].dfGCPPixel &&
pasGCPs[0].dfGCPY == pasGCPs[1].dfGCPY &&
pasGCPs[2].dfGCPY == pasGCPs[3].dfGCPY &&
pasGCPs[0].dfGCPX == pasGCPs[3].dfGCPX &&
pasGCPs[1].dfGCPX == pasGCPs[2].dfGCPX &&
pasGCPs[0].dfGCPY != pasGCPs[2].dfGCPY &&
pasGCPs[0].dfGCPX != pasGCPs[1].dfGCPX
) {
padfGeoTransform[1] =
(pasGCPs[1].dfGCPX - pasGCPs[0].dfGCPX) /
(pasGCPs[1].dfGCPPixel - pasGCPs[0].dfGCPPixel)
padfGeoTransform[2] = 0.0
padfGeoTransform[4] = 0.0
padfGeoTransform[5] =
(pasGCPs[2].dfGCPY - pasGCPs[1].dfGCPY) /
(pasGCPs[2].dfGCPLine - pasGCPs[1].dfGCPLine)
padfGeoTransform[0] =
pasGCPs[0].dfGCPX - pasGCPs[0].dfGCPPixel * padfGeoTransform[1]
padfGeoTransform[3] =
pasGCPs[0].dfGCPY - pasGCPs[0].dfGCPLine * padfGeoTransform[5]
return true
}
/* -------------------------------------------------------------------- */
/* Compute source and destination ranges so we can normalize */
/* the values to make the least squares computation more stable. */
/* -------------------------------------------------------------------- */
let min_pixel = pasGCPs[0].dfGCPPixel
let max_pixel = pasGCPs[0].dfGCPPixel
let min_line = pasGCPs[0].dfGCPLine
let max_line = pasGCPs[0].dfGCPLine
let min_geox = pasGCPs[0].dfGCPX
let max_geox = pasGCPs[0].dfGCPX
let min_geoy = pasGCPs[0].dfGCPY
let max_geoy = pasGCPs[0].dfGCPY
for (let i = 1; i < nGCPCount; ++i) {
min_pixel = Math.min(min_pixel, pasGCPs[i].dfGCPPixel)
max_pixel = Math.max(max_pixel, pasGCPs[i].dfGCPPixel)
min_line = Math.min(min_line, pasGCPs[i].dfGCPLine)
max_line = Math.max(max_line, pasGCPs[i].dfGCPLine)
min_geox = Math.min(min_geox, pasGCPs[i].dfGCPX)
max_geox = Math.max(max_geox, pasGCPs[i].dfGCPX)
min_geoy = Math.min(min_geoy, pasGCPs[i].dfGCPY)
max_geoy = Math.max(max_geoy, pasGCPs[i].dfGCPY)
}
const EPS = 1.0e-12
if (
Math.abs(max_pixel - min_pixel) < EPS ||
Math.abs(max_line - min_line) < EPS ||
Math.abs(max_geox - min_geox) < EPS ||
Math.abs(max_geoy - min_geoy) < EPS
) {
return false // degenerate in at least one dimension.
}
const pl_normalize = new Float64Array(6),
geo_normalize = new Float64Array(6)
pl_normalize[0] = -min_pixel / (max_pixel - min_pixel)
pl_normalize[1] = 1.0 / (max_pixel - min_pixel)
pl_normalize[2] = 0.0
pl_normalize[3] = -min_line / (max_line - min_line)
pl_normalize[4] = 0.0
pl_normalize[5] = 1.0 / (max_line - min_line)
geo_normalize[0] = -min_geox / (max_geox - min_geox)
geo_normalize[1] = 1.0 / (max_geox - min_geox)
geo_normalize[2] = 0.0
geo_normalize[3] = -min_geoy / (max_geoy - min_geoy)
geo_normalize[4] = 0.0
geo_normalize[5] = 1.0 / (max_geoy - min_geoy)
/* -------------------------------------------------------------------- */
/* In the general case, do a least squares error approximation by */
/* solving the equation Sum[(A - B*x + C*y - Lon)^2] = minimum */
/* -------------------------------------------------------------------- */
let sum_x = 0.0
let sum_y = 0.0
let sum_xy = 0.0
let sum_xx = 0.0
let sum_yy = 0.0
let sum_Lon = 0.0
let sum_Lonx = 0.0
let sum_Lony = 0.0
let sum_Lat = 0.0
let sum_Latx = 0.0
let sum_Laty = 0.0
for (let i = 0; i < nGCPCount; ++i) {
const { pdfGeoX: pixel, pdfGeoY: line } = GDALApplyGeoTransform(
pl_normalize,
pasGCPs[i].dfGCPPixel,
pasGCPs[i].dfGCPLine
)
const { pdfGeoX: geox, pdfGeoY: geoy } = GDALApplyGeoTransform(
geo_normalize,
pasGCPs[i].dfGCPX,
pasGCPs[i].dfGCPY
)
sum_x += pixel
sum_y += line
sum_xy += pixel * line
sum_xx += pixel * pixel
sum_yy += line * line
sum_Lon += geox
sum_Lonx += geox * pixel
sum_Lony += geox * line
sum_Lat += geoy
sum_Latx += geoy * pixel
sum_Laty += geoy * line
}
const divisor =
nGCPCount * (sum_xx * sum_yy - sum_xy * sum_xy) +
2 * sum_x * sum_y * sum_xy -
sum_y * sum_y * sum_xx -
sum_x * sum_x * sum_yy
/* -------------------------------------------------------------------- */
/* If the divisor is zero, there is no valid solution. */
/* -------------------------------------------------------------------- */
if (divisor === 0.0) return false
/* -------------------------------------------------------------------- */
/* Compute top/left origin. */
/* -------------------------------------------------------------------- */
const gt_normalized = new Float64Array(6)
gt_normalized[0] =
(sum_Lon * (sum_xx * sum_yy - sum_xy * sum_xy) +
sum_Lonx * (sum_y * sum_xy - sum_x * sum_yy) +
sum_Lony * (sum_x * sum_xy - sum_y * sum_xx)) /
divisor
gt_normalized[3] =
(sum_Lat * (sum_xx * sum_yy - sum_xy * sum_xy) +
sum_Latx * (sum_y * sum_xy - sum_x * sum_yy) +
sum_Laty * (sum_x * sum_xy - sum_y * sum_xx)) /
divisor
/* -------------------------------------------------------------------- */
/* Compute X related coefficients. */
/* -------------------------------------------------------------------- */
gt_normalized[1] =
(sum_Lon * (sum_y * sum_xy - sum_x * sum_yy) +
sum_Lonx * (nGCPCount * sum_yy - sum_y * sum_y) +
sum_Lony * (sum_x * sum_y - sum_xy * nGCPCount)) /
divisor
gt_normalized[2] =
(sum_Lon * (sum_x * sum_xy - sum_y * sum_xx) +
sum_Lonx * (sum_x * sum_y - nGCPCount * sum_xy) +
sum_Lony * (nGCPCount * sum_xx - sum_x * sum_x)) /
divisor
/* -------------------------------------------------------------------- */
/* Compute Y related coefficients. */
/* -------------------------------------------------------------------- */
gt_normalized[4] =
(sum_Lat * (sum_y * sum_xy - sum_x * sum_yy) +
sum_Latx * (nGCPCount * sum_yy - sum_y * sum_y) +
sum_Laty * (sum_x * sum_y - sum_xy * nGCPCount)) /
divisor
gt_normalized[5] =
(sum_Lat * (sum_x * sum_xy - sum_y * sum_xx) +
sum_Latx * (sum_x * sum_y - nGCPCount * sum_xy) +
sum_Laty * (nGCPCount * sum_xx - sum_x * sum_x)) /
divisor
/* -------------------------------------------------------------------- */
/* Compose the resulting transformation with the normalization */
/* geotransformations. */
/* -------------------------------------------------------------------- */
const gt1p2 = new Float64Array(6)
const inv_geo_normalize = new Float64Array(6)
if (!GDALInvGeoTransform(geo_normalize, inv_geo_normalize)) return false
GDALComposeGeoTransforms(pl_normalize, gt_normalized, gt1p2)
GDALComposeGeoTransforms(gt1p2, inv_geo_normalize, padfGeoTransform)
/* -------------------------------------------------------------------- */
/* Now check if any of the input points fit this poorly. */
/* -------------------------------------------------------------------- */
if (bApproxOK) {
// FIXME? Not sure if it is the more accurate way of computing
// pixel size
const dfPixelSize =
0.5 *
(Math.abs(padfGeoTransform[1]) +
Math.abs(padfGeoTransform[2]) +
Math.abs(padfGeoTransform[4]) +
Math.abs(padfGeoTransform[5]))
if (dfPixelSize === 0.0) {
Log.info('dfPixelSize = 0')
return false
}
for (let i = 0; i < nGCPCount; i++) {
const dfErrorX =
pasGCPs[i].dfGCPPixel * padfGeoTransform[1] +
pasGCPs[i].dfGCPLine * padfGeoTransform[2] +
padfGeoTransform[0] -
pasGCPs[i].dfGCPX
const dfErrorY =
pasGCPs[i].dfGCPPixel * padfGeoTransform[4] +
pasGCPs[i].dfGCPLine * padfGeoTransform[5] +
padfGeoTransform[3] -
pasGCPs[i].dfGCPY
if (
Math.abs(dfErrorX) > dfPixelThreshold * dfPixelSize ||
Math.abs(dfErrorY) > dfPixelThreshold * dfPixelSize
) {
Log.info(
`dfErrorX/dfPixelSize = ${
Math.abs(dfErrorX) / dfPixelSize
}, dfErrorY/dfPixelSize = ${Math.abs(dfErrorY) / dfPixelSize}`
)
return false
}
}
}
return true
}
