4DModeller provides some simple map plotting functions that help you create maps from many different types of data.
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.Polygons
We’ll first start off by plotting some polygon data on a map. To do this we’ll use some of the data we use in our Hydrology tutorial, where we display a catchment area on map. First we need to retrieve the example data, then load in the raster data from file.
fdmr::retrieve_tutorial_data(dataset = "hydro")##
## Tutorial data extracted to /home/runner/fdmr/tutorial_data/hydro
norway_polygon_location <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson")Next we need to manipulate the data to ensure we plot the area correctly,
norway_polygon <- rgdal::readOGR(norway_polygon_location)## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /home/runner/.local/share/proj:/usr/share/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.1-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## OGR data source with driver: GeoJSON
## Source: "/home/runner/fdmr/tutorial_data/hydro/Kvilldal_Catch_Boundary.geojson", layer: "Kvilldal_Catch_Boundary"
## with 1 features
## It has 5 fields## Warning in rgdal::readOGR(norway_polygon_location): Z-dimension discarded
norway_polygon <- sf::st_as_sf(norway_polygon,
coords = c("longitude", "latitude"),
crs = "+proj=utm +zone=32"
)
sfc <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84")Now we can plot the polygon on a leaflet map using
plot_map.
fdmr::plot_map(polygon_data = sfc)Rasters
We can also plot raster data onto the map. We’ll now read some raster
data from a NetCDF file. We
do this using the raster::stack function which creates a
RasterStack object which itself contains a number of
RasterLayer objects. These hold ERA5
precipitation data for each day of the month for October 2021.
era5_data_filepath <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "era5_land_daily.nc")
era5_precip <- raster::stack(era5_data_filepath)## Loading required namespace: ncdf4
era5_precip## class : RasterStack
## dimensions : 10, 14, 140, 31 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1000002 (x, y)
## extent : 5.9, 7.3, 59.2, 60.2 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2021.10.01, X2021.10.02, X2021.10.03, X2021.10.04, X2021.10.05, X2021.10.06, X2021.10.07, X2021.10.08, X2021.10.09, X2021.10.10, X2021.10.11, X2021.10.12, X2021.10.13, X2021.10.14, X2021.10.15, ...We can see that we’ve got a stack of 31 layers, one for each day of October. To list all the dates in the stack we can do
names(era5_precip)## [1] "X2021.10.01" "X2021.10.02" "X2021.10.03" "X2021.10.04" "X2021.10.05"
## [6] "X2021.10.06" "X2021.10.07" "X2021.10.08" "X2021.10.09" "X2021.10.10"
## [11] "X2021.10.11" "X2021.10.12" "X2021.10.13" "X2021.10.14" "X2021.10.15"
## [16] "X2021.10.16" "X2021.10.17" "X2021.10.18" "X2021.10.19" "X2021.10.20"
## [21] "X2021.10.21" "X2021.10.22" "X2021.10.23" "X2021.10.24" "X2021.10.25"
## [26] "X2021.10.26" "X2021.10.27" "X2021.10.28" "X2021.10.29" "X2021.10.30"
## [31] "X2021.10.31"As we want to plot this on a Leaflet map we need to ensure the coordinate reference system (CRS) is correct.
projected_raster <- raster::projectRaster(era5_precip, crs = "+proj=utm +zone=32")As we only want to plot the data that overlaps with our rain
catchment area polygon we plotted earlier, we mask the data using raster::mask
and then crop it using terra::crop.
crop_era5 <- raster::mask(projected_raster, norway_polygon)
crop_era5 <- terra::crop(crop_era5, raster::extent(norway_polygon), snap = "near")We’re now read to plot the raster data. Let’s select a single layer from the stack, the first day of October, and plot the precipitation measurements for that day over the polygon we previously plotted.
raster_image <- crop_era5$X2021.10.01Now we’re ready to plot both the polygon and raster data.
fdmr::plot_map(raster_data = raster_image, polygon_data = sfc, palette = "viridis")Interactive mapping
You can plot data on an interactive leaflet map
using our plot_interactive_map function. This launches a Shiny app allowing you to customise
your plot to better investigate your data. Below are a couple of
examples showing how to use the interactive plotting tool.
Hydroelectric Power
We want to create a plot on the map covering the rain catchment area around the power plant. To do this we read in a GeoJSON file containing the coordinates of the catchment area and make sure it’s in the correct coordinate reference system to plot on our map.
norway_polygon_location <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "Kvilldal_Catch_Boundary.geojson")
norway_polygon <- rgdal::readOGR(norway_polygon_location)## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## OGR data source with driver: GeoJSON
## Source: "/home/runner/fdmr/tutorial_data/hydro/Kvilldal_Catch_Boundary.geojson", layer: "Kvilldal_Catch_Boundary"
## with 1 features
## It has 5 fields## Warning in rgdal::readOGR(norway_polygon_location): Z-dimension discarded
norway_polygon <- sf::st_as_sf(norway_polygon,
coords = c("longitude", "latitude"),
crs = "+proj=utm +zone=32"
)
norway_polygon <- sf::st_transform(norway_polygon, crs = "+proj=longlat +datum=WGS84")Next we read ECMWF Reanalysis v5 (ERA5) precipitation data from a NetCDF file.
era5_data_filepath <- fdmr::get_tutorial_datapath(dataset = "hydro", filename = "era5_land_daily.nc")
era5_precip <- raster::stack(era5_data_filepath)
era5_precip <- raster::projectRaster(era5_precip, crs = "+proj=utm +zone=32")The data is now ready to pass to our interactive plotting function.
fdmr::plot_interactive_map(raster_data = era5_precip, polygon_data = norway_polygon)