Gridding and Contour Mapping in R
My blog entry on gridding and contouring.
Notes on how to make contour maps in R
As a retired Geophysicist, I spent a career making contour maps. I have found it to be challenging to make good contour maps in R, and so as part of my own learning process, I have documented the necessary steps in hopes that this may help others involved in the same struggles.
Ultimately I want to create filled contours that represent a surface, based on some random collection of input points, and display those contours on top of a detailed basemap, such as OpenStreetMap.
I will try to use what I believe are the most recent and most likely to survive packages wherever possible. Both static and interactive maps will be generated. I will also be explicit in which packages I am using. I found it rather frustrating to be looking at examples people had on the web, and then wondering which package supplied which function. A rule I read many years ago was “explicit is better than implicit”, and I will try to follow that dictum.
Note that there are several challenges that make this a somewhat difficult exercise.
First of all, this area of R has undergone and continues to see a fair bit of churn. There are a number of packages that address parts of this problem, and they tend to overlap with each other. Googling for help can easily lead to using a superseded package, since there have been so many recent developments. By the same token, the newer packages have not, in all cases, covered the functionality of the older packages, so some use of the older packages may still be necessary.
Secondly, there is a basic issue of coordinate systems. I will display maps in both geographic (that is, Lat-Long) coordinates and in projected (X-Y) coordinates.
Third, interpolating points to a grid is partly art and partly science, so I will give a bit of guidance there.
Fourth, interpolating to a grid in order to then generate contours should be done in a projected X-Y coordinate system, so keeping track of what data is in which system, and projecting back and forth are key.
My test data (seen below) consists of rainfall points in the Houston area from a selection of personal weather stations on August 19, 2022. I will create a simple tibble, and then attach proper geodetic information while converting it to an sf file. I will also transform the sf file to an appropriate projected UTM coordinate system.
Summary of process
- Create geodetically proper data in geographic and/or projected coordinates
- Preliminary work : basemap, Area of Interest
- Plot contours of measurement point density
- Create a grid to define interpolation locations
- Build an interpolator using:
- Inverse weighted nearest neighbor
- Kriging
- Spline fitting
- Thin-plate model
- Interpolate/Predict points onto the grid
- Calculate errors and plot
- Generate contours
- Plot contours over a basemap