The ATPOL grid was created in the late 1960s at the Institute of Botany at Jagiellonian University in Kraków. (Zając 1978a, 1978b) describes the background and methodology. Łukasz Komsta and Marek Verey carried out the extensive mathematical research and GIS implementation (Komsta 2016; Verey 2017). Algorithms provided by Komsta on OpenATPOL are basis for implementation in atpolR package.
In our example we’ll use a distribution map of Erigeron acris L. subsp. acris from (Zając and Zając 2019). The image was scanned at 150 dpi and georeferenced with QGIS.
par(mar = c(0, 0, 0, 0))
<- system.file("extdata/eriacr.tif", package = "atpolR")
tif <- terra::rast(tif)
r ::plotRGB(r) terra
There is hundreds of records. To get them all, we will use
check_atpol_square()
function, which takes the POINT
coordinates and a raster as arguments and checks if the value of raster
cell corresponding to some arbitrary buffer around the POINT equals
zero. As the points are drawn in centers of ATPOL 10 km grid, the values
of atpol10k() centroids with an buffer with radius of 1200 m will be
checked. It may be necessary to adjust a buffer slightly depending of
the quality of scan and precision of georeferencing.
The raster usually consist of 3 layers, one for each R, G, B component. The difference between them are visible on Fig. @ref(fig:rgbraster).
We can see a lot of green and blue components in the scan shown in Fig. @ref(fig:eriacr). We will only look at the first layer for now. Please keep in mind that other layers may be more useful depending on the scan quality and subsequent image processing.
The layer’s values are continuous, ranging from 0 to 255. To simplify
and facilitate analysis, we will classify the layer by assigning only
two values: 0 and 1. We’ll make a reclassification matrix called
rclmat
and use the classify()
function form
the terra package.
<- c(0,120, 0,
m 120,255, 1)
<- matrix(m, ncol=3, byrow=TRUE)
rclmat <- terra::classify(r[[1]], rclmat, include.lowest = TRUE) rc
After preparing the raster, we can load the package and run the
check_atpol_square()
function for all 10k grids:
library(atpolR)
<- atpol10k()|>
eriacr ::mutate(a = mapply(function(x) check_atpol_square(x, rc), centroid)) |>
dplyr::filter(a == "YES") dplyr
Let’s display the results:
eriacr#> Simple feature collection with 709 features and 2 fields
#> Active geometry column: geometry
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 170032.6 ymin: 140353.6 xmax: 839912.3 ymax: 778164.1
#> Projected CRS: ETRS89 / Poland CS92
#> First 10 features:
#> Name geometry centroid a
#> 1 AB15 POLYGON ((220203.9 698774.2... POINT (225196.1 693780.2) YES
#> 2 AB21 POLYGON ((180184 688879.5, ... POINT (185177.4 683882.5) YES
#> 3 AB23 POLYGON ((200188.2 688840.5... POINT (205181.3 683844.9) YES
#> 4 AB73 POLYGON ((200121.5 638968.2... POINT (205117.5 633972.6) YES
#> 5 AB83 POLYGON ((200110.3 628991.5... POINT (205106.9 623996) YES
#> 6 AB94 POLYGON ((210103.7 619001, ... POINT (215100.6 614006.2) YES
#> 7 AB95 POLYGON ((220106.8 618988.3... POINT (225103.5 613994.1) YES
#> 8 AB96 POLYGON ((230109.2 618976, ... POINT (235105.7 613982.3) YES
#> 9 AC05 POLYGON ((220097.7 609011.9... POINT (225094.9 604017.7) YES
#> 10 AC06 POLYGON ((230100.4 609000.3... POINT (235097.4 604006.7) YES
There is 709 observations (grids) in data set.
Let’s have a closer look on BE square:
<- atpol100k() |>
BE subset(Name == "BE") |>
::st_bbox()
sfpar(pty = "s")
plot(NA, type = "n", xlim = c(BE[1], BE[3]), ylim = c(BE[2], BE[4]), axes = FALSE, xlab = "", ylab = "")
::plot(rc, legend = FALSE, add = TRUE)
terra
atpol100k() |>
subset(Name == "BE") |>
::st_cast("LINESTRING") |>
sf::plot(add = TRUE, col = "blue", lwd = 1.2)
terra
|>
eriacr subset(substr(Name, 1, 2) == "BE") |>
::st_set_geometry("centroid") |>
sf::plot(pch = 16, cex = 1.2, col = "blue", add = TRUE) terra
Thus far, so god. Our data set corresponds to published data. In the next step we will add our own observations.
Please keep in mind that check_atpol_square()
function
can produce false positive results in the case of noisy scans. it also
does not recognize the various shapes used in original publications for
various sites.
Assume we want to add our own observations the data set. If we
already know the grid name, we can simply use the Name to filter out the
grid created by the atpol10k()
function. We can use
latlon_to_grid()
function if we don’t have the grid but
have coordinates. In the following example, we are using both
methods:
<- atpol10k() |>
myData ::filter(Name %in% c("BE68",
dplyrlatlon_to_grid(51.13619, 16.95069, 4))) |>
::mutate(a = "myData") dplyr
And let’s add them to above BE square plot:
|>
myData ::st_set_geometry("centroid") |>
sf::plot(pch = 16, cex = 1.8, col = "red", add = TRUE) terra
atpolR package provides a function
plot_points_on_atpol()
which can be used to visualize the
data set. Let’s merge our two data sets (eriacr
and
myData
) and plot them together. For removing any duplicates
we can use unique.data.frame()
function from base R, or
distinct(Name)
from dplyr package.
<- eriacr |>
eriacr rbind(myData) |>
unique.data.frame()
And final plot.
plotPoitsOnAtpol(eriacr$centroid, main = "Erigeron acris subsp. acris", cex = 0.6)
Two basic functions latlon_to_grid()
and
grid_to_latlon()
allows to quickly convert geographical
coordinates (given in WGS 84 latitude and longitude degrees) to ATPOL
grid and from grid to coordinates respectively.
latlon_to_grid(51.01234, 17.23456, 4)
#> [1] "CE50"
latlon_to_grid(51.01234, 17.23456, 6)
#> [1] "CE5086"
The firs two arguments latlon_to_grid()
function are
latitude and longitude respectively. The third argument is the length of
returned grid; it might be even number between 2 and 12.
grid_to_latlon("CE50")
#> [1] 51.04487 17.21736
By default grid_to_latlon()
returns the center of grid
square. If you wish to get it’s corners, you can pass another 2
arguments to the function, which are X and Y offsets, like:
grid_to_latlon("CE50", xoffset = 1, yoffset = 1)
#> [1] 51.00114 17.29032
for bottom right corner.
ATPOL 10km x 10km and 100km x 100k grids are generated by
atpol10k()
and atpol100k()
functions
respectively. It returns set of simple features geometries with grids as
polygons:
atpol100k()
#> Simple feature collection with 44 features and 1 field
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 170009.4 ymin: 110297.1 xmax: 870031.4 ymax: 808651.9
#> Projected CRS: ETRS89 / Poland CS92
#> First 10 features:
#> Name geometry
#> 1 AB POLYGON ((170084 619056.2, ...
#> 2 AC POLYGON ((170014.5 519234.8...
#> 3 AD POLYGON ((170011.8 459334.6...
#> 4 AE POLYGON ((170072 359522.7, ...
#> 5 BA POLYGON ((270217.8 718621.1...
#> 6 BB POLYGON ((270112.1 618931.2...
#> 7 BC POLYGON ((270062.8 519180.9...
#> 8 BD POLYGON ((270060.2 459322.1...
#> 9 BE POLYGON ((270101.2 359577, ...
#> 10 BF POLYGON ((270198.4 259898.5...
For 10k grid it returns a centroids as well:
atpol10k()
#> Simple feature collection with 4400 features and 1 field
#> Active geometry column: geometry
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: 170009.4 ymin: 110297.1 xmax: 870031.4 ymax: 808651.9
#> Projected CRS: ETRS89 / Poland CS92
#> First 10 features:
#> Name geometry centroid
#> 101 AB00 POLYGON ((170215.7 708847.7... POINT (175207.9 703850.1)
#> 102 AB01 POLYGON ((180217.9 708825.6... POINT (185210 703828.7)
#> 103 AB02 POLYGON ((190219.4 708804.1... POINT (195211.3 703808)
#> 104 AB03 POLYGON ((200220 708783.2, ... POINT (205211.9 703787.9)
#> 105 AB04 POLYGON ((210219.9 708763.1... POINT (215211.7 703768.4)
#> 106 AB05 POLYGON ((220219.1 708743.6... POINT (225210.8 703749.6)
#> 107 AB06 POLYGON ((230217.5 708724.8... POINT (235209.1 703731.5)
#> 108 AB07 POLYGON ((240215.3 708706.6... POINT (245206.8 703714)
#> 109 AB08 POLYGON ((250212.4 708689.1... POINT (255203.9 703697.2)
#> 110 AB09 POLYGON ((260208.9 708672.3... POINT (265200.3 703681.1)
Please note, that ATPOL grids are projected in EPSG:2180 coordinate reference system, commonly used in Poland.
Function atpol_div(grid, divider)
divides any given
grid
to smaller grids by 2, 4 or 5 as proposed in (Verey and Komsta 2018).
par(mar = c(0, 0, 0, 0))
<- boundaryPL()
b plot(atpol100k()$geometry)
plot(b, col = "red", add = TRUE)