corny/exdata/openstreetmap.py

@@ -7,6 +7,128 @@ from corny import *
 import pandas
 from shapely.geometry import Point,Polygon
 
+
+class OsmGraph:
+    def __init__(self, network="rivers", bbox=None):
+
+        self.network = "none"
+        if network.startswith("river"):
+            self.network = "rivers"
+        elif network.startswith("road"):
+            self.network = "roads"
+
+        self.bbox = bbox
+        self.graph = None
+        self.shape = None
+
+    def download_graph(self, bbox=None, **kwargs):
+        """
+        Download a network Graph
+        """
+        import osmnx as ox
+        if bbox is None:
+            bbox = self.bbox
+        try:
+            self.graph = ox.graph_from_bbox(
+                bbox[3], bbox[2], bbox[0], bbox[1],
+                **kwargs)
+        except:
+            self.graph = None
+            print("! OsmGraph: could not find a graph for network %s." % self.network)
+            
+        return(self)
+
+    def make_shape(self):
+        """
+        Convert a graph to a shape
+        """
+        if self.graph is None:
+            print("! OsmGraph: missing graph, use download_graph() first.")
+            return()
+        
+        import osmnx as ox
+        try:
+            self.shape = ox.graph_to_gdfs(
+                self.graph,
+                nodes=False, edges=True,
+                fill_edge_geometry=True)
+        except:
+            self.shape = None
+            print("! OsmGraph: cannot convert network graph of %s." % self.network)
+
+        return(self)
+
+    def get(self):
+
+        graph_kw = dict()
+        if self.network == "rivers":
+            graph_kw = dict(
+                simplify         = False,
+                truncate_by_edge = True,
+                retain_all       = True, 
+                custom_filter    = '["waterway"~"river"]'
+            )
+        elif self.network == "roads":
+            graph_kw = dict(
+                simplify         = True,
+                truncate_by_edge = True,
+                retain_all       = True
+            )
+
+        self.download_graph(self.bbox, **graph_kw)
+        self.make_shape()
+        return(self)
+
+    def get_rivers(self, bbox):
+        """
+        Get shape of network graph of rivers
+        """
+        self.download_graph(
+            bbox,
+            simplify         = False,
+            truncate_by_edge = True,
+            retain_all       = True, 
+            custom_filter    = '["waterway"~"river"]'
+            )
+        self.make_shape()
+        return(self.shape)
+
+    def get_roads(self, bbox):
+        """
+        Get shape of network graph of roads
+        """
+        self.download_graph(
+            bbox,
+            simplify         = True,
+            truncate_by_edge = True,
+            retain_all       = True
+            )
+        self.make_shape()
+        return(self.shape)
+
+    def plot(self, ax):
+        """
+        Plotting routine for OsmGraphs. Expects ax.
+        """
+        if self.shape is None:
+            print("! OsmGraph %s: Nothing to plot." % self.network)
+        else:
+            plot_kw = dict()
+            if self.network == "rivers":
+                plot_kw = dict(
+                    zorder=1,
+                    color="blue", alpha=0.1,
+                    lw=2, ls="-")
+            elif self.network == "roads":
+                plot_kw = dict(
+                    zorder=2,
+                    color="k", alpha=0.3,
+                    lw=1, ls=":")
+
+            self.shape.plot(
+                ax=ax, **plot_kw)
+
+
 #######################################################
 ##### Automatic OSM Lookup from Overpass API
 ##### by Erik Nixdorf

corny/figures.py

@@ -26,6 +26,8 @@ from urllib.request import urlopen, Request
 from PIL import Image
 from matplotlib.dates import YearLocator, MonthLocator, DayLocator, WeekdayLocator, HourLocator, MinuteLocator, DateFormatter
 
+from pykrige import OrdinaryKriging
+from pykrige.kriging_tools import write_asc_grid
 
 ##
 
@@ -366,22 +368,31 @@ def add_scalebar(ax, length=1, location=(0.89, 0.04),
     # buffer for scalebar
     buffer = [patheffects.withStroke(linewidth=1, foreground=color)]
     # Plot the scalebar with buffer
-    ax.plot(bar_xs, [sbcy_1, sbcy_1], transform=utm, color='k',
-        linewidth=linewidth, path_effects=buffer)
+    ax.plot(
+        bar_xs, [sbcy_1, sbcy_1],
+        transform=utm,
+        color='k', linewidth=linewidth,
+        path_effects=buffer)
     # buffer for text
     buffer = [patheffects.withStroke(linewidth=1, foreground=color)]
     
     # Plot the scalebar label
-    t0 = ax.text(sbcx_1, sbcy_2, str(length) + ' ' + units, transform=utm,
+    t0 = ax.text(
+        sbcx_1, sbcy_2, str(length) + ' ' + units,
+        transform=utm,
         horizontalalignment='right', verticalalignment='bottom',
         color=color, zorder=2) #path_effects=buffer
     left = x0+(x1-x0)*0.03
     # Plot the N arrow
-    t1 = ax.text(sbcx_1, y_ur, u'\u25B2\nN', transform=utm,
-        horizontalalignment='right', verticalalignment='top',
+    t1 = ax.text(
+        sbcx_1, y_ur, u'\u25B2\nN',
+        transform=utm,
+        horizontalalignment='right',
+        verticalalignment='top',
         color=color, zorder=2)
     # Plot the scalebar without buffer, in case covered by text buffer
-    ax.plot(bar_xs, [sbcy_1, sbcy_1], transform=utm, color=color,
+    ax.plot(bar_xs, [sbcy_1, sbcy_1], transform=utm,
+            color=color,
         linewidth=linewidth, zorder=3)
 
 def guess_ticks_from_lim(a, b, steps=None):
@@ -396,10 +407,10 @@ def guess_ticks_from_lim(a, b, steps=None):
     x_1 = round(float(a)/round_digit)*round_digit
     x_2 = round(float(b)/round_digit)*round_digit
     A = np.arange(x_1, x_2, round_digit)
-    A[(A>=a) & (A<=b)]
+    A = A[(A>=a) & (A<=b)]
     return A
 
-def add_latlon_ticks(ax, steps=0.01):
+def add_latlon_ticks(ax, steps=0.01, grid=True):
     
     # x_1 = round(float(a)/steps)*steps
     # x_2 = round(float(b)/steps)*steps
@@ -411,11 +422,20 @@ def add_latlon_ticks(ax, steps=0.01):
     ys = guess_ticks_from_lim(extent[2], extent[3], steps)
     ax.set_xticks(xs)
     ax.set_yticks(ys)
-    ax.set_xticklabels(["%.3f°" % x for x in xs])
-    ax.set_yticklabels(["%.3f°" % y for y in ys])
+    xlabels = np.array(["%.3f°" % x for x in xs])
+    if len(xlabels) > 7:
+        xlabels[1::3] = ""
+        xlabels[2::3] = ""
+    ylabels = np.array(["%.3f°" % y for y in ys])
+    if len(ylabels) > 7:
+        ylabels[1::3] = ""
+        ylabels[2::3] = ""
+    ax.set_xticklabels(xlabels)
+    ax.set_yticklabels(ylabels)
     ax.set_xlabel(None)#'Longitude (°E)')
     ax.set_ylabel(None)#'Latitude (°N)')
-    ax.grid(color='k', alpha=0.1)
+    if grid:
+        ax.grid(color='k', alpha=0.1)
 
 
 def add_circle(ax, x, y, radius=1,
@@ -437,11 +457,11 @@ def add_colorbar(ax=None, points=None, label=None,
             ):
             
     cb = plt.colorbar(points, ax=ax, **bar_kw)
-    if ticks:
+    if not ticks is None:
         cb.ax.set_yticks(ticks)
-    if ticklabels:
+    if not ticklabels is None:
         cb.ax.set_yticklabels(ticklabels, **ticks_kw)
-    if label:
+    if not label is None:
         cb.set_label(label, **label_kw)
 
 def set_time_ticks(ax=None, how=None, which='major', fmt=None):
@@ -1306,3 +1326,153 @@ class Animation():
     def save(self, filename='test.mp4'):
         self.animation.save(filename=filename, writer="ffmpeg")
 
+
+# OTHER
+
+def make_bbox(x=None, y=None, pad=0):
+
+    if x is None or y is None:
+        default_bbox = [11.030000-pad, 11.06633+pad, 51.640000-pad, 51.670000+pad]
+        print('No x/y given, assuming bbox', default_bbox)
+        return(default_bbox)
+    else:
+        x = list(x)
+        y = list(y)
+        return([min(x)-pad, max(x)+pad, min(y)-pad, max(y)+pad])
+
+
+class Kriging:
+    """
+    Make Kriging interpolation from data[x,y,z] and plot.
+    """
+    def __init__(self, x, y, z,
+                 resolution = 0.01, # 0.001
+                 z_min=0.00, z_max=0.50,
+                 z_cmap="Spectral", z_resolution=0.02):
+
+        self.resolution = resolution
+        # Convert lat/lon resolution into UTM
+        # ...approximation valid in Germany, use proper reprojection elsewhere
+        self.resolution_m = self.resolution/0.001*70 
+
+        self.x = x
+        self.y = y
+        self.z = z
+        self.z_min = z_min
+        self.z_max = z_max
+        self.z_cmap = z_cmap
+        self.z_resolution = z_resolution
+        # show_rivers = True
+        # show_roads  = True
+
+        # Calculate 1D grid for each dimension
+        self.x_grid = np.arange(
+            self.x.min() - self.resolution,
+            self.x.max() + self.resolution*2,
+            self.resolution)
+        self.y_grid = np.arange(
+            self.y.min() - self.resolution,
+            self.y.max() + self.resolution,
+            self.resolution)
+
+        # Calculate mesh grid for countour plot
+        self.x_mesh, self.y_mesh = np.meshgrid(self.x_grid, self.y_grid) 
+
+
+    def execute(self,
+                variogram="spherical" # linear, power, gaussian, spherical, exponential, hole-effect
+        ):
+        
+        krig = OrdinaryKriging(
+            x=self.x, y=self.y, z=self.z,
+            variogram_model = variogram)
+        
+        self.z_grid, self.variance = krig.execute("grid", self.x_grid, self.y_grid)
+
+
+    def plot_variance(self, ax):
+        
+        self.plot(
+            ax, z=self.variance,
+            z_label="Variance of Kriging prediction",
+            contours_kw=dict(
+                cmap="Greys",
+                levels = 20)
+            )
+
+    
+    def plot_values(self, ax):
+        
+        self.plot(
+            ax, z=self.z_grid,
+            z_label="Water content (m³/m³)",
+            contours_kw=dict(
+                cmap="Spectral",
+                levels = np.arange(
+                self.z_min,
+                self.z_max + self.z_resolution,
+                self.z_resolution),
+                vmin = self.z_min,
+                vmax = self.z_max
+            ))
+        
+    plot_z = plot_values
+
+    def plot(self, ax, z=None, 
+             show_z_points=True, z_label="z",
+             contours_kw=dict()):
+        """
+        Usage:
+            with Figure(size=(6,6),
+                projection="flat",
+                extent=bbox,
+                ) as ax:
+
+                Kriging.plot(ax, z)
+        """
+        if z is None:
+            z = self.z_grid
+        if not "cmap" in contours_kw:
+            contours_kw["cmap"] = self.z_cmap
+
+        p_contours = ax.contourf(
+            self.x_mesh, self.y_mesh, z, len(z),
+            **contours_kw
+            )
+
+        # if self.show_rivers and not rivers is None:
+        #     rivers.plot(ax=ax, alpha=0.1, lw=2, color="blue", ls="-",
+        #             zorder=1)
+        # if show_roads and not roads is None:
+        #     roads.plot(ax=ax, alpha=0.3, lw=1, color="k", ls=":",
+        #             zorder=2)
+
+        if show_z_points:    
+            points_kw = dict(
+                c=self.z, s=30, edgecolor="k", lw=0.5,
+                cmap=self.z_cmap, vmin=self.z_min, vmax=self.z_max)
+        else:
+            points_kw = dict(
+                color="k", marker=".", s=1)
+                
+        p_points = ax.scatter(
+            self.x, self.y,
+            **points_kw, zorder=3)
+        
+        add_scalebar(ax, length=1, location=(0.89, 0.04),
+                    color="k", units='km', m_per_unit=1000)
+        add_latlon_ticks(ax, self.resolution, grid=False)
+        
+        add_colorbar(ax, points=p_contours, label=z_label)
+
+        ax.tick_params(top=True, right=True, direction="inout")
+
+
+    def export_asc(self, filename="output.asc"):
+        """
+        Export to ESRI standard ASC grid file
+        """
+        write_asc_grid(
+            self.x_grid, self.y_grid,
+            self.z_grid,
+            filename=filename)
\ No newline at end of file

default.cfg

@@ -601,6 +601,7 @@ make_summary = no
 # Export plots to pdf file (yes/no):
 make_PDF = yes
 # Plot more columns (default is: N, pressure, rel humidity, temperature, abs humidity, moisture)
+# Special plots are: tubes, footprint, diurnal, map_xy, map_tiles, map_interpolation
 PDF_plots = footprint
 # Displayed data range for vol. soil moisture in the plots.
 # Units: m³/m³
@@ -630,6 +631,27 @@ satellite_zoom = 10
 # - watercolor       # Stamen water color
 tiles = osm
 
+# Interpolation options if "map_interpolation" is in "PDF_plots".
+# Variogram model used for Ordinary Kriging interpolation (select one):
+# - spherical   # Spherical model (default)
+# - linear      # Linear model
+# - power       # Power model
+# - gaussian    # Gaussian model
+# - exponential # Exponential model
+variogram = spherical
+# Spatial resolution of the interpolated grid.
+# Units: lat/lon degree, e.g. 0.001 (equals ~70 meters)
+interpolation_ll_resolution = 0.01
+# Resolution of the colorbar.
+# Units: m³/m³
+interpolation_sm_resolution = 0.02
+# Export raster map as ASCII Grid using the ESRI standard (yes/no):
+interpolation_export = yes
+# Add a river network to the plot using OSM (yes/no):
+interpolation_plot_rivers = yes
+# Add a road network to the plot using OSM (yes/no):
+interpolation_plot_roads = yes
+
 ## PNG
 # EXPERIMENTAL: Export a certain plot to a png file and upload it to a webserver directly (yes/no):
 make_PNG = no

examples/rover.cfg

@@ -122,6 +122,10 @@ latitude_column = LatDec
 longitude_column = LongDec
 altitude_column = Alt
 
+## Muons
+# Column for muon counts
+muon_column = 
+
 [clean]
 
 # Cut the data, find outliers, fill gaps, and smooth out noise.
@@ -141,7 +145,11 @@ bbox =
 # Exclude certain areas (e.g. tunnels or buildings). Make sure to add spaces before each line (table):
 # Table: lat, lon, radius (m)
 exclude_table = 
-    
+  
+## UTM conversion
+# Convert lat/lon data to UTM coordinates with EPSG projection. Default is: 31468. Other examples: 25832, ...
+utm_epsg = 31468
+  
 ## Scoring
 # Scoring method (select one):
 # - none # none
@@ -194,7 +202,14 @@ smooth_radius_method = W_r_approx
 # Spatial average over a number of next neighbours (EXPERIMENTAL):
 number_locations = 
 
-## Gridding
+## Regular Gridding
+# Resample all data points to a regular grid based on a UTM coorindate system given by utm_epsg. 
+# Bounding box of the grid. Leave blank for no gridding.
+grid_bbox = 
+# Grid resolution
+grid_resolution = 50
+
+## Custom Gridding
 # Resample all data points to new locations defined in the following master grid.  Leave blank for no spatial resampling
 # Units: csv file
 grid_master_file = 
@@ -269,6 +284,10 @@ temperature_range = -60, 60
 # Units: m³/m³
 sm_range = 0, 0.8
 
+## Muon count range
+# Units: cph
+muon_range = 
+
 
 ## Split tracks
 # Identify single tracks between locations A and B (yes/no)
@@ -332,6 +351,7 @@ new_incoming_column = NM
 # - Schroen et al. (2015)  # Amplitude scaling by gamma
 # - Rotunno and Zreda (2014)  # linear to local cutoff rigidity, gamma ~ Rc
 # - Hawdon et al. (2014)  # Function of local and reference cutoff rigidity, gamma ~ (Rc - Rc_ref)
+# - McJannet and Desilets et al. (2023)  # Rescales JUNG based on lat lon, Rc
 incoming_method = Zreda et al. (2012)
 # Correction parameter to scale the amplitude of the incoming radiation used for Schroen et al. (2015):
 gamma = 1
@@ -361,6 +381,12 @@ NM_station = JUNG
 # - 1day  #  1 day     temporal resolution
 NM_resolution = 1min
 
+## Muon counts
+# Muon correction parameters based on Stevanato et al. (2022)
+muon_beta = 625
+muon_alpha = 0.0021
+muon_T_ref = 0
+
 [conversion]
 
 # New column name for the soil moisture product
@@ -419,6 +445,7 @@ land_use_forest_corr_err = 0.0
 # - Baatz et al. (2015)  # scale neutrons by a percentage of biomass
 # - Jakobi et al. (2018)  # make use of the epithermal/thermal ratio
 veg_method = none
+
 # Local biomass at the site.
 # Units: kg/m²
 biomass = 0
@@ -575,6 +602,7 @@ make_summary = no
 # Export plots to pdf file (yes/no):
 make_PDF = yes
 # Plot more columns (default is: N, pressure, rel humidity, temperature, abs humidity, moisture)
+# Special plots are: tubes, footprint, diurnal, map_xy, map_tiles, map_interpolation
 PDF_plots = footprint_length, footprint, tubes, map_grid, map_tiles, map_interpolation
 # Displayed data range for vol. soil moisture in the plots.
 # Units: m³/m³
@@ -593,7 +621,7 @@ PDF_time_format = H:M
 PDF_line_style = points
 # Satellite picture zoom level. Warning: using high zoom levels for large extent maps could crash the Google map service. Better use low resolution, e.g., zoom = 10, for faster performance.
 # Units: Google zoom level 1--18
-satellite_zoom = 15
+satellite_zoom = 16
 # Map tiles to use for mapping (select one):
 # - osm              # OpenStreetMap (default)
 # - google           # Google Maps
@@ -602,7 +630,28 @@ satellite_zoom = 15
 # - stamen           # Stamen terrain 
 # - toner            # Stamen toner
 # - watercolor       # Stamen water color
-tiles = osm
+tiles = satellite-google
+
+# Interpolation options if "map_interpolation" is in "PDF_plots".
+# Variogram model used for Ordinary Kriging interpolation (select one):
+# - spherical   # Spherical model (default)
+# - linear      # Linear model
+# - power       # Power model
+# - gaussian    # Gaussian model
+# - exponential # Exponential model
+variogram = spherical
+# Spatial resolution of the interpolated grid.
+# Units: lat/lon degree, e.g. 0.001 (equals ~70 meters)
+interpolation_ll_resolution = 0.001
+# Resolution of the colorbar.
+# Units: m³/m³
+interpolation_sm_resolution = 0.02
+# Export raster map as ASCII Grid using the ESRI standard (yes/no):
+interpolation_export = yes
+# Add a river network to the plot using OSM (yes/no):
+interpolation_plot_rivers = yes
+# Add a road network to the plot using OSM (yes/no):
+interpolation_plot_roads = yes
 
 ## PNG
 # EXPERIMENTAL: Export a certain plot to a png file and upload it to a webserver directly (yes/no):

examples/spatial_kriging.py

+# %%
+"""
+This example shows how to apply spatial Kriging interpolation
+and plot neat maps based on data from Corny output.
+"""
+import pandas
+from corny.figures import Figure, make_bbox, Kriging
+from corny.exdata.openstreetmap import OsmGraph
+
+# %%
+"""
+Read data from Corny example output (created using rover.cfg).
+"""
+data = pandas.read_csv(
+    "output/rover-2019082112.csv",
+    index_col=0, parse_dates=True)
+
+# Reduce to the relevant pieces
+data = data[["utmx","utmy","lat","lon","moisture"]].dropna()
+
+# %%
+"""
+Apply Kriging
+"""
+# Setup up Kriging instance
+OKI = Kriging(data["lon"], data["lat"], data["moisture"],
+            resolution=0.001,
+            z_min=0.00, z_max=0.50,
+            z_resolution=0.02)
+
+# Run Kriging operator, creates OKI.z_grid
+OKI.execute()
+
+# Export as ESRI ASCII grid file
+OKI.export_asc("output/sm_kriging_grid.asc")
+
+# %%
+"""
+Make bounding box from data
+"""
+resolution = 0.001
+bbox = make_bbox(
+    data["lon"].values,
+    data["lat"].values,
+    pad=resolution)
+
+# %%
+"""
+Download map features
+"""
+Rivers = OsmGraph("rivers", bbox).get()
+Roads  = OsmGraph("roads",  bbox).get()
+
+# %%
+"""
+Create the plot
+"""
+with Figure(
+    size   = (9,9),
+    layout = (2,1),
+    projection = "flat",
+    extent = bbox,
+    abc=[
+        "Kriging map of water content",
+        "Quality of Kriging"
+        ]
+    ) as axes:
+
+    ax = axes[0]
+    # Plot interpolated soil moisture values
+    OKI.plot_values(ax)
+    # Plot map features
+    Rivers.plot(ax)
+    Roads.plot(ax)
+ 
+    ax = axes[1]
+    # Plot variance map (quality of interpolation)
+    OKI.plot_variance(ax)
+    # Plot map features
+    Rivers.plot(ax)
+    Roads.plot(ax)
+
+# %%

install/python_packages/requirements-conda-pip.txt

@@ -14,4 +14,5 @@ kivymd
 kivy-garden-mapview
 mapview
 kivy-garden
-kivy-garden.graph
\ No newline at end of file
+kivy-garden.graph
+pykrige
\ No newline at end of file

install/python_packages/requirements-linux.txt

@@ -35,4 +35,5 @@ kivy
 kivymd
 mapview
 kivy-garden-mapview
-kivy-garden.graph
\ No newline at end of file
+kivy-garden.graph
+pykrige
\ No newline at end of file

install/python_packages/requirements-windows-2.txt

@@ -30,3 +30,4 @@ kivymd
 mapview
 kivy-garden-mapview
 kivy-garden.graph
+pykrige
\ No newline at end of file

instantPASDy.py

@@ -551,7 +551,7 @@ def main(configfile=None, cfg_prefix=None):
 
         # if no geographic information available we stop processing here
         if data[lat].isna().sum()==len(data) or data[lon].isna().sum()==len(data):
-            print('! Cannot find any GPS data )\':')
+            print('! Cannot find any GPS data:')
             raise ValueError('No GPS data available, processing stopped.')
 
         # replace altitude column with data from DEM
@@ -2614,12 +2614,21 @@ def main(configfile=None, cfg_prefix=None):
                 smrange = xsplit(config['output']['PDF_sm_range'], range=2, type=float)
             else:
                 smrange = (0,1)
-            time_range = (data.index.max()-data.index.min()).total_seconds()/60/60/24/365
+            
+            time_range_hours = (data.index.max()-data.index.min()).total_seconds()/60/60
+            time_range_days  = time_range_hours/24
+            time_range_years = time_range_days/365
+            time_range       = time_range_years
+            
             pdf_plots_list = xsplit(config['output']['PDF_plots'])
 
             # Date range in plots, minimum a full day
-            date1 = datetime(data.index.min().year, data.index.min().month, data.index.min().day, 0,0,0)
-            date2 = datetime(data.index.max().year, data.index.max().month, data.index.max().day, 23,59,59)
+            if time_range_hours > 12:
+                date1 = datetime(data.index.min().year, data.index.min().month, data.index.min().day, 0,0,0)
+                date2 = datetime(data.index.max().year, data.index.max().month, data.index.max().day, 23,59,59)
+            else:
+                date1 = datetime(data.index.min().year, data.index.min().month, data.index.min().day, data.index.min().hour, 0,0)
+                date2 = datetime(data.index.max().year, data.index.max().month, data.index.max().day, data.index.max().hour, 59,59)
 
             # New
             from corny.reports import Report
@@ -2649,8 +2658,10 @@ def main(configfile=None, cfg_prefix=None):
                 kw_fig_time = dict(x_major_ticks="months", x_minor_ticks="days", x_major_fmt="%d\n%b", x_minor_fmt="%d")
             elif time_range > 6/24/365:
                 kw_fig_time = dict(x_major_ticks="days", x_minor_ticks="hours", x_major_fmt="%d\n%b", x_minor_fmt="%H:%M")
-            else:
+            elif time_range > 2/24/365:
                 kw_fig_time = dict(x_major_ticks="hours", x_major_fmt="%H")
+            else:
+                kw_fig_time = dict(x_major_ticks="hours", x_major_fmt="%H:%M", x_minor_ticks="minutes")
             
             with Report(figure_name2, "Data Report: %s" % title) as R:
                 with Progress("| Making PDF report", replace=True) as Pr:
@@ -3556,6 +3567,100 @@ def main(configfile=None, cfg_prefix=None):
                                         ax.scatter( poi.lon, poi.lat, c='black', marker='+', s=100, transform=ccrs.PlateCarree())
 
                             R.add_image(F.buff)
+                    
+
+                    if "map_interpolation" in pdf_plots_list:
+                        
+                        var = sm
+
+                        Pr.update("Map interpolation")
+                        if not data[var].isnull().all() and lat != '' and lat in data.columns:
+
+                            from corny.figures import make_bbox, Kriging
+                            from corny.exdata.openstreetmap import OsmGraph
+
+                            # Reduce data to the relevant pieces
+                            data_x = data[[lat,lon,var]].dropna()
+
+                            # Config
+                            resolution   = gc(config, "output", "interpolation_ll_resolution", alt=0.01, dtype=float)
+                            z_resolution = gc(config, "output", "interpolation_sm_resolution", alt=0.02, dtype=float)
+                            variogram    = gc(config,"output","variogram", alt="spherical")
+
+                            # Setup up Kriging instance
+                            OKI = Kriging(
+                                data_x[lon], data_x[lat],
+                                data_x[var],
+                                resolution = resolution,
+                                z_min=smrange[0], z_max=smrange[1],
+                                z_resolution=z_resolution)
+
+                            # Run Kriging operator, creates OKI.z_grid
+                            OKI.execute(variogram=variogram)
+
+                            # Export as ESRI ASCII grid file
+                            if cc(config, "output", "interpolation_export"):
+                                OKI.export_asc(output_basename + '.asc')
+                                file_save_msg(
+                                    output_basename + '.asc',
+                                    name="ASCII Grid",
+                                    end="\n")
+
+                            # Make bounding box from data
+                            bbox = make_bbox(
+                                data_x[lon].values,
+                                data_x[lat].values,
+                                pad=resolution)
+
+                            # Download map features
+                            if cc(config, "output", "interpolation_plot_rivers"):
+                                Rivers = OsmGraph("rivers", bbox).get()
+                            if cc(config, "output", "interpolation_plot_roads"):
+                                Roads  = OsmGraph("roads",  bbox).get()
+
+                            R.add_page()
+                            R.add_title("Map (interpolated) of %s" % var)
+                            R.add_par("Kriging map of water content created with Ordinary Kriging and a %s variogram model. The resolution of the interpolated grid has been set to %.4f° (~%.1f meters)." \
+                                      % (variogram, OKI.resolution, OKI.resolution_m))
+
+                            # Create the plot
+                            with Figure(
+                                size   = (9,9),
+                                layout = (1,1),
+                                projection = "flat",
+                                extent = bbox,
+                                save="buff", format="png"
+                                ) as F:
+
+                                ax = F.axes
+                                # Plot interpolated soil moisture values
+                                OKI.plot_values(ax)
+                                # Plot map features
+                                Rivers.plot(ax)
+                                Roads.plot(ax)
+                            
+                            R.add_image(F.buff)
+                            R.add_page()
+                            R.add_title("Map (interpolated) of Variance")
+                            R.add_par("Quality of the Kriging interpolation expressed as variance. The variance is lowest near the data points and increases with distance, indicating the increasing uncertainty of the interpolation map with distance from the actual data.") 
+
+                            with Figure(
+                                size   = (9,9),
+                                layout = (1,1),
+                                projection = "flat",
+                                extent = bbox,
+                                save="buff", format="png"
+                                ) as F:
+
+                                ax = F.axes
+                                # Plot variance map (quality of interpolation)
+                                OKI.plot_variance(ax)
+                                # Plot map features
+                                Rivers.plot(ax)
+                                Roads.plot(ax)
+                            
+                            R.add_image(F.buff)
+
 
                     for pdf_plots_var in pdf_plots_list: