[Download] [Project page] [Online map] [Tutorial (Polish)]

*Świetnie opracowany po polsku opis siatki dla nie-programistów można odnaleźć na stronie Instytutu Botaniki PAN.*

This project is a reference, well-tested implementation of conversion algorithms for Polish geobotanical grid called ATPOL.

The grid was developed in Institute of Botany, Jagiellonian University, Krakow, Poland in 1970s for the immortal work called *"Distribution atlas of vascular plants in Poland"*. It covers 700 x 700 km area (with whole Poland inside).

The grid was initially drawn by hand on available maps, then copied also by hand, as no tools like xero were then available. When GPS tools became popular, Institute of Botany published its online version with computed coordinates (with 1 second accuracy).

These coordinates were used by researchers as a reference, however the mathematical background of this grid was never published. In Polish internet fora, there was a long discussion without any definite conclusion. The only one possible conversion was to take the tabularized grid data and to perform the linear (incorrect, but the only one possible) interpolation.

The extensive mathematical research done in 2015 uncovered the mathematical mystery of the ATPOL grid. The exact solution was found by reverse engineering of the tabularized values and fitting dozens of geographical projections. The findings and proposed algorithm was published (in Polish) in:

and this is the main citing reference for OpenATPOL. As one of the reviewers of above paper explicitly asked for consulting these results with Institute of Botany before publication, the author was forced to try and got the confirmation in private correspondence, that this reasoning is correct. The further numerical investigation in press was done by Marek Verey. **Finally, on 16th October 2017, the creators of ATPOL grid confirmed officially the reasoning described in these two papers, confirming also the use of WGS84 coordinate system.**

Although the projection seems to be quite strange (as there were no printed maps in such projection, moreover the assumed Earth radius is equal to 6390 km), this projection was chosen as a compromise by the ATPOL authors to fit all available hand-drawn maps. Therefore, these algorithms can be surely used without any risk that they are incorrect.

The ATPOL projection is simple conic **centrographic** projection with one tangent (the cone touches the sphere but does not cross it) at the latitude \(\varphi_0 = 52^\circ\)N. The central (vertical) meridian is chosen to be \(\lambda_0 = 19^\circ\)E. The central point of the projection \((52^\circ N, 19^\circ E)\) is defined to have coordinates \(x_0 = 330\)km, \(y_0 = 350\)km. The Earth radius is chosen to be \(R = 6390\)km. The coordinates increase towards South and East (they are equal to zero in NW corner of the grid). It is quite similar to Equidistant Conic Projection (EQDC), but the latter is orthographic (the projection rays goes in ATPOL from the Earth center, whereas in EQDC they are orthogonal and parallel).

The first two letters denote the square with 100 km size, as in the figure above. The smaller 10 x 10km square units are denoted by two additional numbers, but in reverse way (\(y\) is the first coordinate). The division can be done recursively, so AB1234 is a square with 1 km dimension, AB987654 is 100m, and so on.

The \(x,y\) coordinates are never used in practice, however they represent an exact position in ATPOL projection and they are an intermediate step when geographical coordinates are converted to ATPOL grid. When grid is given, the position inside the square can be done in relative way by giving \(o_x\) and \(o_y\) values in range \((0,1)\). Two zeros means upper left (NW) corner of the square, whereas two 0.5 values means the middle of the square.

The coordinates are in currently used WGS84 coordinate system. Although the old maps could be drawn in Pulkowo 1942 coordinate systems, the difference between these coordinate systems is less than hand-copying error of the maps that time.

For example, the coordinates of the Palace of Culture and Science in Warsaw can be represented as:

- As "lonlat" pair, \(\varphi = 21.006062^\circ\)N, \(\lambda = 52.231727^\circ\)E.
- As "xy" pair, \(x = 467.0110005\)km, \(y = 322.2659527\)km.
- As a point inside ED (100 x 100 km) square, with relative position inside \(o_x = 0.67011001, o_y = 0.22265953\).
- As a point inside ED26 (10 x 10 km) square, with relative position inside \(o_x = 0.7011001, o_y = 0.2265953\).
- As a point inside ED2627 (100 x 100 km) square, with relative position inside \(o_x = 0.011001, o_y = 0.2265953\).
- As a point inside ED262720 (100 x 100 m) square, with relative position inside \(o_x = 0.11001, o_y =0.65953\).
- As a point inside ED26272061 (10 x 10 m) square, with relative position inside \(o_x = 0.1001, o_y = 0.5953\).
- As a point inside ED2627206161 (1 x 1 m) square, with relative position inside \(o_x = 0.001, o_y = 0.953\).
- (and so on).

One should remember about the GPS accuracy, so the last example has any sense only with geodesic GPS receivers.

The \(o_x\) and \(o_y\) values are never given in publications, they are only needed to store exact position and give the ability to perform the reverse computations.

**There is apparently no simple solution** to set centrographic tangent conic projection in GIS software with PROJ4 library.

The most similar solution can be set with the following custom projection:

```
+proj=eqdc +lat_1=52 +lat_2=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000
+b=6390000 +ellps=sphere +x_0=330000 +y_0=-350000
```

or, when "+axis=esu" is not supported:

```
+proj=eqdc +lat_1=52 +lat_2=52 +lat_0=52 +lon_0=19 +a=6390000
+b=6390000 +ellps=sphere +x_0=330000 +y_0=350000
```

The second case will give reversed \(y\) coordinates, counting from lower left (SW) corner.

*However, this is only an approximation and results in \(\approx\) 400 m error on the grid edges.*

VBA implementation (for Excel and OpenOffice) is available for download. Other languages passed the tests and the snippets can be browsed in SVN tree.

For each language, two tests for accuracy check are designed and available:

- The short test, which computes the points given as a reference below and compares with the reference values. The sum of errors is returned, it should not exceed 1e-10.
- The long test, which generates pseudorandom (repeatable) 10 000 pairs of "xy" points and converts each point to latitudes and longitudes. For each grid length in range 2,4,...,10,12 each point is converted to grid and back to coordinates 10 times. Maximal error is reported (it should be around 1e-15). This test is extensive and allowed to catch several floating point rounding errors and deal with them. Beware that the computation is very long, the spreadsheet
**can hang**so do not run this function before backup or besides any production data.

Void functions, which require to pass pointers to target variables to store the results. The grid char pointer must point to an allocated place to store grid results there.

```
void atpol_latlon_to_xy(double lat, double lon, double *x, double *y);
void atpol_xy_to_latlon(double x, double y, double *lat, double *lon);
void atpol_xy_to_grid(double x, double y, char *grid, int len, double *xoffset, double *yoffset);
void atpol_grid_to_xy(const char *grid, double xoffset, double yoffset, double *x, double *y);
void atpol_grid_to_latlon(const char *grid, double xoffset, double yoffset, double *lat, double *lon);
void atpol_latlon_to_grid(double lat, double lon, char *grid, int length, double *xoffset, double *yoffset);
```

This implementation tries to avoid OOP as much as possible to force consistency with other languages. Methods are static and return object arrays, containing the results.

```
public static Object[]
latlon_to_xy(double lat, double lon) // 0 - x (double), 1 - y (double)
public static Object[]
xy_to_latlon(double x, double y) // 0 - lat (double), 1 - lon (double)
public static Object[]
xy_to_grid(double x, double y, int len) // 0 - grid (string), 1 - xoffset (double), 2 - yoffset (double)
public static Object[]
grid_to_xy(String grid, double xoffset, double yoffset) // 0 - x (double), 1 - y (double)
public static Object[]
latlon_to_grid(double lat, double lon, int len) // 0 - grid (string), 1 - xoffset (double), 2 - yoffset (double)
public static Object[]
grid_to_latlon(String grid, double xoffset, double yoffset) // 0 - lat (double), 1 - lon (double)
```

```
var atpol = {
latlon_to_xy: function (lat, lon) // returns [x,y]
xy_to_latlon: function (x, y) // returns [lat,lon]
xy_to_grid: function (x, y, len) // returns [grid, xoffset, yoffset]
grid_to_xy: function (grid, xoffset, yoffset) // returns [x,y]
grid_to_latlon: function (grid, xoffset, yoffset) // returns [lat,lon]
latlon_to_grid: function (lat, lon, length) // returns [grid, xoffset, yoffset]
}
```

```
class atpol {
function latlon_to_xy($lat, $lon) // returns array($x, $y)
function xy_to_latlon($x, $y) // returns array($lat, $lon);
function xy_to_grid($x, $y, $len=8) // returns array($grid, $xoffset, $yoffset);
function grid_to_xy($grid, $xoffset, $yoffset) // returns array($x, $y)
function grid_to_latlon($grid, $xoffset, $yoffset) // returns array($lat, $lon);
function latlon_to_grid($lat, $lon, $length) // returns array($grid, $xoffset, $yoffset);
}
```

(tested with freepascal compiler)

```
unit atpol;
interface
procedure atpol_latlon_to_xy(lat : double; lon: double; var x : double; var y : double);
procedure atpol_xy_to_latlon(x: double; y: double; var lat: double; var lon: double);
procedure atpol_xy_to_grid(x : double; y: double; length : integer; var grid: string;
var xoffset: double; var yoffset: double);
procedure atpol_grid_to_xy(grid : string; xoffset : double; yoffset : double;
var x : double; var y : double);
procedure atpol_latlon_to_grid(lat : double; lon : double; var grid : string;
length : integer; var xoffset : double; var yoffset : double);
procedure atpol_grid_to_latlon(grid : string; xoffset : double; yoffset : double;
var lat : double; var lon : double);
```

```
latlon_to_xy($lat, $lon) # return($xx,$yy);
xy_to_latlon($x, $y) # return($lat,$lon);
xy_to_grid($x,$y,$len) # return($grid, $xoffset, $yoffset);
grid_to_xy($grid, $xoffset, $yoffset) # return($x,$y);
latlon_to_grid($lat, $lon, $len) # return($lat,$lon);
grid_to_latlon($grid, $xoffset, $yoffset) # return($lat,$lon);
```

```
def latlon_to_xy(lat, lon) # return (x, y)
def xy_to_latlon(x, y) # return (lat, lon)
def xy_to_grid(x, y, len=8) # return (grid, xoffset, yoffset)
def grid_to_xy(grid, xoffset, yoffset) # return (x, y)
def grid_to_latlon(grid, xoffset, yoffset) # return (lat, lon)
def latlon_to_grid(lat, lon, length) # return (grid, xoffset, yoffset)
```

VBA implementation allows to perform the coordinate conversions inside the spreadsheets:

**ATPOL_GRID_TO_LATLON**(grid, xoffset, yoffset)

-returns an array (lat, lon). The offset parameters are optional, default values are 0.5 (the middle of the square).**ATPOL_GRID_TO_XY**(grid, xoffset, yoffset)

-returns an array (x, y). The offset parameters are optional, default values are 0.5 (the middle of the square).**ATPOL_LATLON_TO_GRID**(lat, lon, length)

-returns an array (grid, xoffset, yoffset, x, y). The length parameter is optional (default = 8, i.e. 100 meters square).**ATPOL_LATLON_TO_XY**(lat, lon)

-returns an array (x, y).**ATPOL_XY_TO_GRID**(x, y, length)

-returns an array (grid, xoffset, yoffset, x, y). The length parameter is optional (default = 8, i.e. 100 meters square).**ATPOL_XY_TO_LATLON**(x, y)

-returns an array (lat, lon)**ATPOL_TEST**()

-runs the short test, returns sum of errors**ATPOL_TEST_LONG**()

-runs the long test, returns maximal error**(can hang the spreadsheet program, do not run before data backup)**.

The algorithm for "latlon" pair to "xy" pair conversion:

- \(\lambda = \lambda - \lambda_0\)
- \(r = \frac{\cos \varphi_0}{\sin \varphi_0} - \tan (\varphi - \varphi_0)\)
- \(x = r \sin (\lambda \sin\varphi_0)\)
- \(y = r \sin (\lambda \cos\varphi_0)\)
- \(y = y - \frac{\cos \varphi_0}{\sin \varphi_0}\)
- \(x = Rx + x_0\)
- \(y = Ry + x_0\)

The algorithm for "xy" pair to "latlon" pair conversion:

- \(x = \frac{x - x_0}{R}\)
- \(y = \frac{y - y_0}{R}\)
- \(y = y + \frac{\cos \varphi_0}{\sin \varphi_0}\)
- \(\lambda = \frac{\arctan(x/y)}{\sin \varphi_0}\)
- \(r = \sqrt{x^2+y^2}\)
- \(\varphi = \varphi_0 - \arctan\left(r-\frac{\cos \varphi_0}{\sin \varphi_0}\right)\)
- \(\lambda = \lambda + \lambda_0\)

The next two algorithms seem to be crazy, however this is the way to avoid jumping of the results between several squares due to floating point error during division by any power of 10 (yes, it happened very rarely, but tests caught this problem).

The algorithm for "xy" pair to "grid" conversion:

- Convert \(x\) and \(y\) to meters and round it to create an integer values \(x_s\) and \(y_s\)
- Convert to 6-character string with optional leading zeros.
- Create a grid string, by taking the first character of \(x_s\) and the first of \(y_s\), increasing the ASCII value by 17 (zero becomes A and so on).
- For \(i = 1 \dots 5\) append to the grid \(i\)-th character of \(y_s\) followed by \(i\)-th character of \(x_s\).
- Create \(o_x\) and \(o_y\) strings, containing the \(x\) and \(y\) values expressed as milimeters (multiply by 1e6) and convert to string with length equal to 9, optionally panned with leading zeros.
- Take the last \(9 - l/2\) characters of \(o_x\) and \(o_y\) and add zero with decimal point at the beginning.
- Convert \(o_x\) and \(o_y\) from string to float.
- Return \(l\) first characters of grid together with \(o_x\) and \(o_y\).

The algorithm for "grid" to "xy" conversion:

- Create \(x_s\) by taking the first character of grid, decreasing ASCII value by 17 (A becomes 0 and so on).
- Create \(y_s\) by taking the second character of grid, decreasing ASCII value by 17 (A becomes 0 and so on).
- For \(i = 1 \dots l/2\): append \(x_s\) with the \((2i+1)\)-th character of the grid, then append \(y_s\) with the \(2i\)-th character of the grid.
- Append zeros to \(x_s\) and \(y_s\) until their length is equal to 6.
- Convert \(x_s\) and \(y_s\) from string to int.
- \(x = x_s/1000 + 10^{3-l/2}o_x\)
- \(y = y_s/1000 + 10^{3-l/2}o_y\)

These reference values (beware of exponents after "b" letter) were computed up to 90 significant digits using Maxima Computer Algebra System. They were also checked in R using MPFR multiple precision library interface.

```
lat = 55, lon = 24
x = 6.500315410 9413219363 4638921909
2919102347 6549426346 4294406114
1134263122 2586709523 6877638936 b2
y = 4.106161777 0643609028 0069398239
9868109797 7416098516 8475411589
7741848959 9526507236 9127814945 b0
lat = 49, lon = 15
x = 3.707418900 7307473069 6146050746
2762635775 1822453853 9452623925
1805565421 1002353757 3371334723 b1
y = 6.768262355 9270039774 4316813206
3136538237 1028407494 1210768315
2116602614 4742493968 5399610358 b2
lat = 49, lon = 24
x = 6.960533606 1617843913 9749832254
1755018117 4703675369 6342716067
6503001532 8499992312 2712744907 b2
y = 6.722945679 5827199940 2311874332
7574763977 2116443227 5129393298
2056341816 4582670512 9387468615 b2
lat = 52, lon = 19
x = 330
y = 350
x = 0, y = 0
lat = 5.503040399 3648806391 6131639786
9569363580 3219644367 1837742899
2831279682 1573747622 5554628753 b1
lon = 1.384022731 8521004431 6191689836
7708417395 0606591263 3069285468
7864255005 5905203650 9154275744 b1
x = 700, y = 0
lat = 5.500351550 5218481835 1436251940
2356517141 5525549010 2366181839
2910549810 8234368660 1539449084 b1
lon = 2.478270718 4271129766 1229102707
7864236604 5958438506 2611201880
3643885511 3577052279 0628302380 b1
x = 0, y = 700
lat = 4.877384783 4747808675 2723438685
7484285450 8438794922 2721477719
9175666700 1434651537 3381662738 b1
lon = 1.451445359 4615022781 0015146712
5841166221 7940635007 0428788960
0558335581 1832098353 1992814714 b1
x = 700, y = 700
lat = 4.875047607 0495021286 9950531623
1417070334 8694846514 1777654287
4672439505 7853796537 0474833601 b1
lon = 2.402761076 3560529927 7468614203
6788418600 5147161668 5964092529
8571408752 0740760263 6494787753 b1
```

OpenATPOL is created and maintaned by Łukasz Komsta.

Released under GPLv3 license. Last update: 2017.11.12