Analysis
The analysis area is the Krajkowo nature reserve located in the Greater Poland voivodeship. It was established in 1958 in order to protect the breeding places of birds, especially the grey heron and the great black cormorant, and to protect the landscape of the Warta oxbow.
Vector data
Data on nature reserves can be found in General Geographic Databases.
We can obtain them using the geodb_download() function.
Let’s do that.
If you run into any problem with the download, remember that you can
pass another download method from download.file() as a
function argument.
The downloaded database consists of many files in the GML
(Geography Markup Language) format. A brief description of the
structure of this database can be found here.
The table with the nature reserves is in the
“PL.PZGIK.201.30__OT_TCRZ_A.xml” file. We can use the
sf package and its read_sf() function to
load it.
Let’s check the structure of our data.
In simple terms, it is a spatial table consisting of 110 observations (rows) and 28 variables (columns). The names of the objects are located in the nazwa column, which allow us to select the Krajkowo reserve only.
We can display it in two basic ways:
- Using the
plot()function and directly specifying the column with the object geometry:plot(krajkowo$geometry) - Using the
plot()andst_geometry()functions that obtain geometry from the vector layer. In the first case, we need to know the name of the column with geometries (e.g.geometry,geom, etc.), while in the second case, the geometry is selected automatically (it is the safer and preferable way).
We can also calculate the area of this polygon.
## 165.2744 [ha]The st_area() function returned the area in m^2, and
after the conversion we got the result of 165 ha.
Raster data
Now let’s move on to the stage of downloading the orthophotomap. We
use the ortho_request() function that show us which images
are available for the analyzed area. We need to provide our Krajkowo
polygon as the argument of this function.
We can display the resulting table using the code below.
## sheetID year resolution composition sensor CRS
## 1 N-33-142-B-d-4-2 2004 0.50 B/W Analog PL-1992
## 2 N-33-142-B-d-4-4 2004 0.50 B/W Analog PL-1992
## 3 N-33-142-B-d-4-4 2010 0.25 RGB Digital PL-1992
## 4 N-33-142-B-d-4-2 2010 0.25 RGB Digital PL-1992
## 5 N-33-142-B-d-4-4 2010 0.25 CIR Digital PL-1992
## 6 N-33-142-B-d-4-2 2010 0.25 CIR Digital PL-1992
## 7 N-33-142-B-d-4-2 2016 0.25 CIR Digital PL-1992
## 8 N-33-142-B-d-4-2 2016 0.25 RGB Digital PL-1992
## 9 N-33-142-B-d-4-4 2016 0.25 CIR Digital PL-1992
## 10 N-33-142-B-d-4-4 2016 0.25 RGB Digital PL-1992To complete our task, we need to obtain a near infrared data. So in
the next step, we select those rows for which the
composition column has the value of “CIR”.
Then let’s sort the table according to the year the photo was taken, with the most recent images at the top of the table.
Let’s display the table again and select the newest compositions.
## sheetID year resolution composition sensor seriesID
## 7 N-33-142-B-d-4-2 2016 0.25 CIR Digital 69837
## 9 N-33-142-B-d-4-4 2016 0.25 CIR Digital 69837
## 22 N-33-142-B-d-4-2 2013 0.25 CIR Digital 69903
## 24 N-33-142-B-d-4-4 2013 0.25 CIR Digital 69903
## 5 N-33-142-B-d-4-4 2010 0.25 CIR Digital 69763
## 6 N-33-142-B-d-4-2 2010 0.25 CIR Digital 69763Note that the result has a pair of objects (images). This means that
our Krajkowo reserve is depicted in two photos within one series.
Therefore, the seriesID column is used to combine smaller
images into a larger mosaic.
## sheetID year resolution composition sensor seriesID
## 7 N-33-142-B-d-4-2 2016 0.25 CIR Digital 69837
## 9 N-33-142-B-d-4-4 2016 0.25 CIR Digital 69837The tile_download() function is used to download
orthophotomaps by taking our selected table as the main argument. We can
also specify the output folder with the outdir
argument.
## 1/2
## 2/2Processing
Let’s load the downloaded orthophotomaps using the
read_stars() function from the stars
package, which allows working with spatiotemporal arrays. In our case,
it is a raster consisting of three bands (NIR, R, G) for only one point
in time. We don’t need to load the entire data array into memory - we
can read the file’s metadata instead by using the proxy
argument.
img1 = read_stars("data/69837_329609_N-33-142-B-d-4-2.TIF", proxy = TRUE)
img2 = read_stars("data/69837_329613_N-33-142-B-d-4-4.TIF", proxy = TRUE)Now we can perform two operations: rasters merging and cropping to
the reserve area. The use of a proxy allows to get the
result almost immediately, while processing the entire image
(proxy = FALSE) would take several minutes. Images have
their own specific Coordinate Reference Systems, so let’s make sure it
is the correct one after merging. It should be EPSG 2180 in this
case.
img = st_mosaic(img1, img2)
st_crs(img) = 2180 # overwrite CRS to be sure
img = st_crop(img, krajkowo)Let’s display the effect using the plot() function, and
define the input bands with the rgb argument. It creates a
composition consisting of three bands: NIR, R and G in our case. The
composition below is shown in infrared, not in natural colors, which may
be misinterpreted from the rgb argument name.
In the last step, we calculate the NDVI using the near infrared (1) and red (2) bands.
calc_ndvi = function(img) (img[1] - img[2]) / (img[1] + img[2])
ndvi = st_apply(img, MARGIN = c("x", "y"), FUN = calc_ndvi)
plot(ndvi, main = "NDVI", col = hcl.colors(10, palette = "RdYlGn"))
A surprising observation is the relatively low NDVI values for the forest area. There are two reasons for this, i.e. the photos are taken in mid-March (before the start of the growing season) and probably have not been calibrated. For this reason, a better source of data for analysis may be satellite images, which are calibrated spectrally and obtained continuously (if no cloudiness occurs).