In the article I will tell you about the results of my "research", compose an algorithm for converting an arbitrary raster map into tiles that are understandable for applications and, along the way, get acquainted with such concepts as an ellipsoid, datum, coordinate system, projection.
Our Earth does not have the shape of a ball, and not even the shape of an ellipsoid, this complex figure is called a geoid. The fact is that the masses inside the Earth are not evenly distributed, so in some places the Earth is slightly concave, in others it is slightly convex. If we take the territory of a separate country or continent, then its surface with sufficient accuracy can be aligned with the surface of some ellipsoid, if the center of this ellipsoid is slightly shifted along three coordinate axes relative to the Earth's center of mass. Such an ellipsoid is called a reference ellipsoid, it is suitable for describing only the local area of the Earth for which it was created. At large distances from this site, it can have a very large discrepancy with the Earth's surface. An ellipsoid whose center coincides with the center of mass of the Earth is called a common terrestrial ellipsoid. Clear,that the reference ellipsoid better describes its local part of the Earth than the general terrestrial, but the general terrestrial is suitable for the entire surface of the Earth.
To describe the ellipsoid, only two independent values are sufficient: the equatorial radius (usually denoted by a) and the polar radius (b), but instead of the second independent value, polar contraction f = (ab) / a is usually used. This is the first thing we need in our algorithm as an object - an ellipsoid. For different parts of the Earth in different years, different researchers have calculated many reference ellipsoids, information about them is given in the form of data: a (in meters) and 1 / f (dimensionless). Oddly enough, for the common terrestrial ellipsoid, there are also many different variants (different a, f), but the difference is not very strong, it is mainly due to the difference in the methods for determining a and f.
struct Ellipsoid {
char *name;
double a; /* () */
double b; /* () */
double al; /* (a-b)/a */
double e2; /* (a^2-b^2)/a^2 */
};
struct Ellipsoid Ellipsoid_WGS84 = {
.name = "WGS84",
.a = 6378137.0,
.al = 1.0 / 298.257223563,
};
struct Ellipsoid Ellipsoid_Krasovsky = {
.name = "Krasovsky",
.a = 6378245.0,
.al = 1.0 / 298.3,
};
The example shows two ellipsoids: the common terrestrial WGS84, used in satellite navigation, and the Krasovsky reference ellipsoid, applicable to Europe and Asia.
Consider another interesting point, the fact is that the shape of the Earth is slow, but changing, so the ellipsoid, which today describes the surface well, in a hundred years may be far from reality. This has little to do with the common terrestrial ellipsoid, since deviations within its same error, but relevant for the reference ellipsoid. Here we come to another concept - a datum. Datum is a set of parameters of an ellipsoid (a, f), its displacements inside the Earth (for a reference ellipsoid), fixed at a certain point in time. More precisely, the datum may not necessarily describe an ellipsoid, it can be parameters of a more complex figure, for example, a quasigeoid.
Surely the question has already arisen: how to move from one ellipsoid or datum to another? For this, each ellipsoid must have a system of geodetic coordinates: latitude and longitude (phi, lambda), the transition is carried out by translating coordinates from one coordinate system to another.
There are various formulas for transforming coordinates. You can first translate geodesic coordinates in one coordinate system into three-dimensional coordinates X, Y, Z, perform shifts and rotations with them, and then convert the resulting three-dimensional coordinates into geodetic coordinates in another coordinate system. You can do it directly. Because all formulas are infinite converging series, then they are usually limited to a few members of the series to achieve the required accuracy. As an example, we will use the Helmert transforms, these are transformations using a transition to three-dimensional coordinates, they consist of the three stages described above. For transformations, in addition to two ellipsoids, you need 7 more parameters: three shifts along three axes, three rotation angles and a scale factor. As you might guess, all parameters can be extracted from datums.But in the algorithm we will not use such an object as a datum, but instead we will introduce a transition object from one coordinate system to another, which will contain: references to two ellipsoids and 7 transformation parameters. This will be the second object of our algorithm.
struct HelmertParam {
char *src, *dst;
struct Ellipsoid *esp;
struct Ellipsoid *edp;
struct {
double dx, dy, dz;
double wx, wy, wz;
double ms;
} p;
//
double a, da;
double e2, de2;
double de2__2, dxe2__2;
double n, n__2e2;
double wx_2e2__ro, wy_2e2__ro;
double wx_n__ro, wy_n__ro;
double wz__ro, ms_e2;
};
struct HelmertParam Helmert_SK42_WGS84 = {
.src = "SK42",
.dst = "WGS84",
.esp = &Ellipsoid_Krasovsky,
.edp = &Ellipsoid_WGS84,
// SK42->PZ90->WGS84 ( 51794-2001)
.p = {23.92, -141.27, -80.9, 0, -0.35, -0.82, -0.12e-6},
};
The example shows the parameters for converting from the Pulkovo 1942 coordinate system to the WGS84 coordinate system. The transformation parameters themselves are a separate topic, for some coordinate systems they are open, for others they are selected empirically, therefore their values may slightly differ in different sources.
In addition to parameters, a transformation function is also needed, it can be direct and for transformation in the opposite direction, we only need a transformation in the opposite direction. I will skip tons of math and give my optimized function.
void setupHelmert(struct HelmertParam *pp) {
pp->a = (pp->edp->a + pp->esp->a) / 2;
pp->da = pp->edp->a - pp->esp->a;
pp->e2 = (pp->edp->e2 + pp->esp->e2) / 2;
pp->de2 = pp->edp->e2 - pp->esp->e2;
pp->de2__2 = pp->de2 / 2;
pp->dxe2__2 = pp->de2__2 + pp->e2 * pp->da / pp->a;
pp->n = 1 - pp->e2;
pp->n__2e2 = pp->n / pp->e2 / 2;
pp->wx_2e2__ro = pp->p.wx * pp->e2 * 2 * rad(0,0,1);
pp->wy_2e2__ro = pp->p.wy * pp->e2 * 2 * rad(0,0,1);
pp->wx_n__ro = pp->p.wx * pp->n * rad(0,0,1);
pp->wy_n__ro = pp->p.wy * pp->n * rad(0,0,1);
pp->wz__ro = pp->p.wz * rad(0,0,1);
pp->ms_e2 = pp->p.ms * pp->e2;
}
void translateHelmertInv(struct HelmertParam *pp,
double lat, double lon, double h, double *latp, double *lonp) {
double sin_lat, cos_lat;
double sin_lon, cos_lon;
double q, n;
if (unlikely(!pp)) {
*latp = lat;
*lonp = lon;
return;
}
sin_lat = sin(lat);
cos_lat = cos(lat);
sin_lon = sin(lon);
cos_lon = cos(lon);
q = 1 / (1 - pp->e2 * sin_lat * sin_lat);
n = pp->a * sqrt(q);
*latp = lat
- ((n * (q * pp->de2__2 + pp->dxe2__2) * sin_lat + pp->p.dz) * cos_lat
- (pp->p.dx * cos_lon + pp->p.dy * sin_lon) * sin_lat
) / (n * q * pp->n + h)
+ (pp->wx_2e2__ro * sin_lon - pp->wy_2e2__ro * cos_lon)
* (cos_lat * cos_lat + pp->n__2e2)
+ pp->ms_e2 * sin_lat * cos_lat;
*lonp = lon
+ ((pp->p.dx * sin_lon - pp->p.dy * cos_lon) / (n + h)
- (pp->wx_n__ro * cos_lon + pp->wy_n__ro * sin_lon) * sin_lat
) / cos_lat
+ pp->wz__ro;
}
Where does all this come from? In a more understandable language, formulas can be found in the proj4 project, but since I needed optimization for the speed of execution, then any calculations of the sine of the sum of the angles were transformed by the formulas, the exponentiations were optimized by veneings in brackets, and the constants were calculated separately.
Now, to get closer to completing the original task of creating tiles, we need to consider a coordinate system called WebMercator. This coordinate system is used in the OsmAnd application and in the web, for example, in Google maps and in OpenStreetMap. WebMercator is a Mercator projection built on a sphere. The coordinates in this projection are the coordinates of the pixel X, Y and they depend on the Z scale, for a zero scale the entire earth's surface (up to about 85 degrees of latitude) is placed on one 256x256 pixel tile, the X, Y coordinates change from 0 to 255, starting from the left top corner, for scale 1 - already 4 tiles, X, Y - from 0 to 511 and so on.
The following functions are used to convert from WebMercator to WGS84:
void XYZ_WGS84(unsigned x, unsigned y, int z, double *latp, double *lonp) {
double s = M_PI / ((1UL << 7) << z);
*lonp = s * x - M_PI;
*latp = asin(tanh(M_PI - s * y));
}
void WGS84_XYZ(int z, double lat, double lon, unsigned *xp, unsigned *yp) {
double s = ((1UL << 7) << z) / M_PI;
*xp = uint_round((lon + M_PI) * s);
*yp = uint_round((M_PI - atanh(sin(lat))) * s);
}
And at the end of the first part of the article, we can already sketch out the beginning of our algorithm for creating a tile. Each tile of 256x256 pixels is addressed by three values: x, y, z, the relationship with coordinates X, Y, Z is very simple: x = (int) (X / 256); y = (int) (Y / 256); z = Z;
void renderTile(int z, unsigned long x, unsigned long y) {
int i, j;
double wlat, wlon;
double lat, lon;
for (i = 0; i < 255; ++i) {
for (j = 0; j < 255; ++j) {
XYZ_WGS84(x * 256 + j, y * 256 + i, z, &wlat, &wlon);
translateHelmertInv(&Helmert_SK42_WGS84, wlat, wlon, 0, &lat, &lon);
/* lat,lon - 42 */
}
}
}
Coordinates in SK42 are already transformed coordinates into the map coordinate system, now it remains to find a pixel on the map by these coordinates and enter its color into a tile pixel at coordinates j, i. This will be the next article, in which we will talk about geodesic projections and affine transformations.