Geo-referenced basemaps

In simulariatools 3.1.0, we enhanced the wrapper function contourPlot2() to interpret and utilize geographical information from a geo-referenced basemap file.

The primary motivation for this new feature is to minimize errors that may arise from potential misalignment between the basemap and the data matrix to be plotted. In previous versions, it was possible for the contour plot to be overlaid on the wrong basemap. Now, by using a GeoTIFF basemap, we can safely assume that the contour plot will be accurately aligned with it without providing the plot bounding box via xlim and ylim arguments. Nonetheless, we must still be careful with the Coordinate Reference System (CRS) of the data. Since the data.frame() to be plotted is not CRS aware, we need to be sure that the coordinates x and y are in the same reference system as the GeoTIFF file.

Here is a simple example of this new feature in action. Let’s consider a sample dataset of yearly average odour units of an imaginary source, produced by the dispersion model CALPUFF. To import the data into the project, we use importSurferGrd() from simulariatools even though there are other ways to accomplish this task. The computational domain of the numerical simulation has lower left coordinates (388000, 4969000) in the UTM32 CRS, an extent of 6000 m in both directions and a grid size of 200 m.

calpuff_file <- "./ou_8760hr_testcase.grd"
ou_testcase <- importSurferGrd(calpuff_file, k = 1000)
#> [1] Z min =  0.3609
#> [1] Z max =  4.697
head(ou_testcase)
#>        x       y       z
#> 1 388000 4969000 0.38817
#> 2 388200 4969000 0.43803
#> 3 388400 4969000 0.48059
#> 4 388600 4969000 0.51562
#> 5 388800 4969000 0.54536
#> 6 389000 4969000 0.57245

There are many ways to obtain a geo-referenced basemap file. Here we download the image from the Italian geoportal by using the downloadBasemap() function of simulariatools. In its most recent version, the function automatically tries to download a geo-referenced TIFF file (GeoTIFF).

downloadBasemap(
    file = "./basemap.tiff",
    xSW = 388000,
    ySW = 4969000,
    xExt = 6000,
    yExt = 6000
)

In order to analyze the file content, it is convenient to use the rast() function of the terra package:

terra::rast("basemap.tiff")
#> class       : SpatRaster
#> size        : 1024, 1024, 3  (nrow, ncol, nlyr)
#> resolution  : 5.859375, 5.859375  (x, y)
#> extent      : 388000, 394000, 4969000, 4975000  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632)
#> source      : basemap.tiff
#> colors rgb  : 1, 2, 3
#> names       : basemap_1, basemap_2, basemap_3

As we can see, this basemap exactly overlaps with the dispersion simulation computational domain.

Now we can quickly plot the odour levels using the bounding box defined by the GeoTIFF file:

contourPlot2(
    ou_testcase,
    basemap = "basemap.tiff",
    nticks = 7,
    levels = c(1, 2, 3, 4, 5),
    fill = FALSE
)

In this case, the bounding box defined by the basemap file and the extent defined by the ou_testcase data.frame are the same. There is no need to explicitly specify the plot bounding box, nor to assume that the basemap has the correct extent.

We can use the new xlim and ylim arguments to crop the contour plot:

contourPlot2(
    ou_testcase,
    basemap = "basemap.tiff",
    xlim = c(389000, 393000),
    ylim = c(4970000, 4974000),
    levels = c(1, 2, 3, 4, 5),
    fill = FALSE
)

As we can easily check by looking at the odour level contour lines, the plot is still correctly aligned with the underlying image, without the need to pass through a new basemap file.

To test the new function even further, we can use a basemap covering a larger geographical area. In this case, we use an already prepared GeoTIFF file created by QGIS with a map from OpenStreetMap, with an extent of 40 km × 40 km:

terra::rast("basemap_large.tiff")
#> class       : SpatRaster
#> size        : 2481, 2481, 4  (nrow, ncol, nlyr)
#> resolution  : 16.12253, 16.12253  (x, y)
#> extent      : 370000, 410000, 4960000, 5000000  (xmin, xmax, ymin, ymax)
#> coord. ref. : WGS 84 / UTM zone 32N (EPSG:32632)
#> source      : basemap_large.tiff
#> names       : basemap_large_1, basemap_large_2, basemap_large_3, basemap_large_4

contourPlot2() reads the extent of the larger basemap and correctly plots the data in ou_testcase on the map:

contourPlot2(
    ou_testcase,
    basemap = "basemap_large.tiff",
    levels = c(1, 2, 3, 4, 5)
) +
    labs(x = NULL, y = NULL)