Skip to content

Commit bde296c

Browse files
kylebarronkylebarron
authored
feat: Extract TileMatrixSet definition from GeoTIFF object (#50)
Port @vincentsarago 's PR here developmentseed/morecantile#187 to JS for use with geotiff.js. This creates our subset of the TMS spec from a COG image. Note that we diverge from the morecantile PR on ordering: we reverse our ordering of `tileMatrices` to be from coarsest to finest.
1 parent 4d388df commit bde296c

File tree

10 files changed

+1251
-7
lines changed

10 files changed

+1251
-7
lines changed

‎packages/deck.gl-cog/src/cog-tile-matrix-set.ts‎

Lines changed: 216 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,216 @@
1+
import type {
2+
TileMatrix,
3+
TileMatrixSet,
4+
} from "@developmentseed/deck.gl-raster";
5+
import type { GeoTIFF, GeoTIFFImage } from "geotiff";
6+
import {
7+
extractGeotransform,
8+
getGeoTIFFProjection,
9+
} from "./geotiff-reprojection";
10+
import proj4, { ProjectionDefinition } from "proj4";
11+
import Ellipsoid from "./ellipsoids.js";
12+
import type { PROJJSONDefinition } from "proj4/dist/lib/core";
13+
14+
// 0.28 mm per pixel
15+
// https://docs.ogc.org/is/17-083r4/17-083r4.html#toc15
16+
const SCREEN_PIXEL_SIZE = 0.00028;
17+
18+
/**
19+
*
20+
* Ported from Vincent's work here:
21+
* https://github.com/developmentseed/morecantile/pull/187/changes#diff-402eedddfa30af554d03750c8a82a09962b44b044976c321b774b484b98e8f48R665
22+
*
23+
* @return {TileMatrixSet}[return description]
24+
*/
25+
export async function parseCOGTileMatrixSet(
26+
tiff: GeoTIFF,
27+
): Promise<TileMatrixSet> {
28+
const fullResImage = await tiff.getImage();
29+
30+
if (!fullResImage.isTiled) {
31+
throw new Error("COG TileMatrixSet requires a tiled GeoTIFF");
32+
}
33+
34+
const imageCount = await tiff.getImageCount();
35+
const bbox = fullResImage.getBoundingBox();
36+
const fullImageWidth = fullResImage.getWidth();
37+
const fullImageHeight = fullResImage.getHeight();
38+
39+
const crs = await getGeoTIFFProjection(fullResImage);
40+
41+
if (crs === null) {
42+
throw new Error(
43+
"Could not determine coordinate reference system from GeoTIFF geo keys",
44+
);
45+
}
46+
47+
const parsedCrs = parseCrs(crs);
48+
49+
const projectToWgs84 = proj4(crs, "EPSG:4326").forward;
50+
const projectTo3857 = proj4(crs, "EPSG:3857").forward;
51+
52+
const boundingBox: TileMatrixSet["boundingBox"] = {
53+
lowerLeft: [bbox[0]!, bbox[1]!],
54+
upperRight: [bbox[2]!, bbox[3]!],
55+
};
56+
57+
const transform = extractGeotransform(fullResImage);
58+
59+
if (transform[1] !== 0 || transform[3] !== 0) {
60+
// TileMatrixSet assumes orthogonal axes
61+
throw new Error(
62+
"COG TileMatrixSet with rotation/skewed geotransform is not supported",
63+
);
64+
}
65+
66+
const cellSize = Math.abs(transform[0]);
67+
68+
const tileWidth = fullResImage.getTileWidth();
69+
const tileHeight = fullResImage.getTileHeight();
70+
71+
const tileMatrices: TileMatrix[] = [
72+
{
73+
// Set as highest resolution / finest level
74+
id: String(imageCount - 1),
75+
scaleDenominator:
76+
(cellSize * metersPerUnit(parsedCrs)) / SCREEN_PIXEL_SIZE,
77+
cellSize,
78+
pointOfOrigin: [transform[2], transform[5]],
79+
tileWidth: fullResImage.getTileWidth(),
80+
tileHeight: fullResImage.getTileHeight(),
81+
matrixWidth: Math.ceil(fullImageWidth / tileWidth),
82+
matrixHeight: Math.ceil(fullImageHeight / tileHeight),
83+
geotransform: transform,
84+
},
85+
];
86+
87+
// Starting from 1 to skip full res image
88+
for (let imageIdx = 1; imageIdx < imageCount; imageIdx++) {
89+
const image = await tiff.getImage(imageIdx);
90+
91+
if (!image.isTiled) {
92+
throw new Error("COG TileMatrixSet requires a tiled GeoTIFF");
93+
}
94+
95+
const tileMatrix = createOverviewTileMatrix({
96+
id: String(imageCount - 1 - imageIdx),
97+
image,
98+
fullWidth: fullImageWidth,
99+
fullHeight: fullImageHeight,
100+
baseTransform: transform,
101+
parsedCrs,
102+
});
103+
tileMatrices.push(tileMatrix);
104+
}
105+
106+
// Reverse to have coarsest level first
107+
tileMatrices.reverse();
108+
109+
return {
110+
crs,
111+
boundingBox,
112+
tileMatrices,
113+
projectToWgs84,
114+
projectTo3857,
115+
};
116+
}
117+
118+
/**
119+
* Coefficient to convert the coordinate reference system (CRS)
120+
* units into meters (metersPerUnit).
121+
*
122+
* From note g in http://docs.opengeospatial.org/is/17-083r2/17-083r2.html#table_2:
123+
*
124+
* > If the CRS uses meters as units of measure for the horizontal dimensions,
125+
* > then metersPerUnit=1; if it has degrees, then metersPerUnit=2pa/360
126+
* > (a is the Earth maximum radius of the ellipsoid).
127+
*/
128+
// https://github.com/developmentseed/morecantile/blob/7c95a11c491303700d6e33e9c1607f2719584dec/morecantile/utils.py#L67-L90
129+
function metersPerUnit(parsedCrs: ProjectionDefinition): number {
130+
switch (parsedCrs.units) {
131+
case "metre":
132+
case "meter":
133+
case "meters":
134+
return 1;
135+
case "foot":
136+
return 0.3048;
137+
case "US survey foot":
138+
return 1200 / 3937;
139+
}
140+
141+
if (parsedCrs.units === "degree") {
142+
// 2 * π * ellipsoid semi-major-axis / 360
143+
const { a } = Ellipsoid[parsedCrs.ellps as keyof typeof Ellipsoid];
144+
return (2 * Math.PI * a) / 360;
145+
}
146+
147+
throw new Error(`Unsupported CRS units: ${parsedCrs.units}`);
148+
}
149+
150+
/**
151+
* Create tile matrix for COG overview
152+
*/
153+
function createOverviewTileMatrix({
154+
id,
155+
image,
156+
fullWidth,
157+
baseTransform,
158+
parsedCrs,
159+
}: {
160+
id: string;
161+
image: GeoTIFFImage;
162+
fullWidth: number;
163+
fullHeight: number;
164+
baseTransform: [number, number, number, number, number, number];
165+
parsedCrs: ProjectionDefinition;
166+
}): TileMatrix {
167+
const width = image.getWidth();
168+
const height = image.getHeight();
169+
170+
// For now, just use scaleX
171+
// https://github.com/developmentseed/morecantile/pull/187/changes#r2621314673
172+
const scaleX = fullWidth / width;
173+
174+
// const scaleY = fullHeight / height;
175+
// if (Math.abs(scaleX - scaleY) > 1e-3) {
176+
// throw new Error("Non-uniform overview scaling detected (X/Y differ)");
177+
// }
178+
179+
const scale = scaleX;
180+
181+
const geotransform: [number, number, number, number, number, number] = [
182+
baseTransform[0] * scale,
183+
baseTransform[1] * scale,
184+
baseTransform[2], // x origin stays the same
185+
baseTransform[3] * scale,
186+
baseTransform[4] * scale,
187+
baseTransform[5], // y origin stays the same
188+
];
189+
const cellSize = Math.abs(geotransform[0]);
190+
191+
const tileWidth = image.getTileWidth();
192+
const tileHeight = image.getTileHeight();
193+
194+
return {
195+
id,
196+
scaleDenominator: (cellSize * metersPerUnit(parsedCrs)) / SCREEN_PIXEL_SIZE,
197+
cellSize,
198+
pointOfOrigin: [geotransform[2], geotransform[5]],
199+
tileWidth,
200+
tileHeight,
201+
matrixWidth: Math.ceil(width / tileWidth),
202+
matrixHeight: Math.ceil(height / tileHeight),
203+
geotransform,
204+
};
205+
}
206+
207+
function parseCrs(crs: PROJJSONDefinition): ProjectionDefinition {
208+
// If you pass proj4.defs a projjson, it doesn't parse it; it just returns the
209+
// input.
210+
//
211+
// Instead, you need to assign it to an alias and then retrieve it.
212+
213+
const key = "__deck.gl-cog-internal__";
214+
proj4.defs(key, crs);
215+
return proj4.defs(key);
216+
}

‎packages/deck.gl-cog/src/ellipsoids.ts‎

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
/**
2+
* Vendored and edited from proj4.js.
3+
*
4+
* In the implementation of metersPerUnit while generated a COG TileMatrixSet,
5+
* we need to know the size of the semi major axis in the case the CRS is in
6+
* degrees.
7+
*
8+
* https://github.com/proj4js/proj4js/blob/e90e5fa6872a1ffc40edb161cbeb4bd5e3bd9db5/lib/constants/Ellipsoid.js
9+
*/
10+
11+
const ellipsoids = {
12+
MERIT: {
13+
a: 6378137,
14+
},
15+
SGS85: {
16+
a: 6378136,
17+
},
18+
GRS80: {
19+
a: 6378137,
20+
},
21+
IAU76: {
22+
a: 6378140,
23+
},
24+
airy: {
25+
a: 6377563.396,
26+
b: 6356256.91,
27+
},
28+
APL4: {
29+
a: 6378137,
30+
},
31+
NWL9D: {
32+
a: 6378145,
33+
},
34+
mod_airy: {
35+
a: 6377340.189,
36+
b: 6356034.446,
37+
},
38+
andrae: {
39+
a: 6377104.43,
40+
},
41+
aust_SA: {
42+
a: 6378160,
43+
},
44+
GRS67: {
45+
a: 6378160,
46+
},
47+
bessel: {
48+
a: 6377397.155,
49+
},
50+
bess_nam: {
51+
a: 6377483.865,
52+
},
53+
clrk66: {
54+
a: 6378206.4,
55+
b: 6356583.8,
56+
},
57+
clrk80: {
58+
a: 6378249.145,
59+
},
60+
clrk80ign: {
61+
a: 6378249.2,
62+
b: 6356515,
63+
},
64+
clrk58: {
65+
a: 6378293.645208759,
66+
},
67+
CPM: {
68+
a: 6375738.7,
69+
},
70+
delmbr: {
71+
a: 6376428,
72+
},
73+
engelis: {
74+
a: 6378136.05,
75+
},
76+
evrst30: {
77+
a: 6377276.345,
78+
},
79+
evrst48: {
80+
a: 6377304.063,
81+
},
82+
evrst56: {
83+
a: 6377301.243,
84+
},
85+
evrst69: {
86+
a: 6377295.664,
87+
},
88+
evrstSS: {
89+
a: 6377298.556,
90+
},
91+
fschr60: {
92+
a: 6378166,
93+
},
94+
fschr60m: {
95+
a: 6378155,
96+
},
97+
fschr68: {
98+
a: 6378150,
99+
},
100+
helmert: {
101+
a: 6378200,
102+
},
103+
hough: {
104+
a: 6378270,
105+
},
106+
intl: {
107+
a: 6378388,
108+
},
109+
kaula: {
110+
a: 6378163,
111+
},
112+
lerch: {
113+
a: 6378139,
114+
},
115+
mprts: {
116+
a: 6397300,
117+
},
118+
new_intl: {
119+
a: 6378157.5,
120+
b: 6356772.2,
121+
},
122+
plessis: {
123+
a: 6376523,
124+
},
125+
krass: {
126+
a: 6378245,
127+
},
128+
SEasia: {
129+
a: 6378155,
130+
b: 6356773.3205,
131+
},
132+
walbeck: {
133+
a: 6376896,
134+
b: 6355834.8467,
135+
},
136+
WGS60: {
137+
a: 6378165,
138+
},
139+
WGS66: {
140+
a: 6378145,
141+
},
142+
WGS7: {
143+
a: 6378135,
144+
},
145+
WGS84: {
146+
a: 6378137,
147+
},
148+
sphere: {
149+
a: 6370997,
150+
b: 6370997,
151+
},
152+
};
153+
154+
export default ellipsoids;

0 commit comments

Comments
 (0)