corny/exdata/openstreetmap.py

@@ -3,12 +3,36 @@ CoRNy exdata.OpenStreetMap
     Get Data from OpenStreetMap
 """
 
+import numpy as np
 from corny import *
 import pandas
 from shapely.geometry import Point,Polygon
 
 
 class OsmGraph:
+
+    road_colors = dict(
+        footway = "C7",
+        mixed = "black",
+        path = "C2",
+        primary = "C3",
+        service = "C8",
+        track = "C5",
+        trunk = "C4",
+        trunk_link = "C6",
+        unclassified = "#CCCCCC",
+        motorway = "red",
+        motorway_link = "pink",
+        residential = "C9",
+        secondary = "C1",
+        tertiary = "C0",
+        road = "#888888",
+        primary_link = "violet",
+        secondary_link = "plum",
+        tertiary_link = "thistle",
+        none = "white"
+        )
+    
     def __init__(self, network="rivers", bbox=None):
 
         self.network = "none"
@@ -19,7 +43,8 @@ class OsmGraph:
 
         self.bbox = bbox
         self.graph = None
-        self.shape = None
+        self.data = None
+        self.road_colors = OsmGraph.road_colors
 
     def download_graph(self, bbox=None, **kwargs):
         """
@@ -38,79 +63,92 @@ class OsmGraph:
             
         return(self)
 
-    def make_shape(self):
+    def to_data(self):
         """
-        Convert a graph to a shape
+        Convert a graph to a data
         """
         if self.graph is None:
-            print("! OsmGraph: missing graph, use download_graph() first.")
-            return()
+            self.data = None
+            # print("! OsmGraph: missing graph, use download_graph() first.")
+            return(self)
         
         import osmnx as ox
         try:
-            self.shape = ox.graph_to_gdfs(
+            self.data = 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)
+            self.data = None
+            # print("! OsmGraph: cannot convert network graph of %s." % self.network)
 
         return(self)
 
     def get(self):
-
-        graph_kw = dict()
+        """
+        Wrapper to generate river or road data
+        """
         if self.network == "rivers":
-            graph_kw = dict(
-                simplify         = False,
-                truncate_by_edge = True,
-                retain_all       = True, 
-                custom_filter    = '["waterway"~"river"]'
-            )
+            self.get_rivers()
+
         elif self.network == "roads":
-            graph_kw = dict(
-                simplify         = True,
-                truncate_by_edge = True,
-                retain_all       = True
-            )
+            self.get_roads()
 
-        self.download_graph(self.bbox, **graph_kw)
-        self.make_shape()
         return(self)
 
-    def get_rivers(self, bbox):
+    def get_rivers(self):
         """
-        Get shape of network graph of rivers
+        Get data of network graph of rivers
         """
         self.download_graph(
-            bbox,
+            self.bbox,
             simplify         = False,
             truncate_by_edge = True,
             retain_all       = True, 
             custom_filter    = '["waterway"~"river"]'
             )
-        self.make_shape()
-        return(self.shape)
+        self.to_data()
+
+        return(self)
 
-    def get_roads(self, bbox):
+    def get_roads(self):
         """
-        Get shape of network graph of roads
+        Get data of network graph of roads
         """
         self.download_graph(
-            bbox,
+            self.bbox,
             simplify         = True,
             truncate_by_edge = True,
             retain_all       = True
             )
-        self.make_shape()
-        return(self.shape)
+        self.to_data()
+
+        if not self.data is None:
+            self.data["highway_str"] = ""
+            for i, row in self.data.iterrows():
+                row = row.copy()
+                if isinstance(row["highway"], list):
+                    if 'primary' in row["highway"]:
+                        self.data.loc[i, 'highway_str'] = 'primary'
+                    elif 'track' in row["highway"]:
+                        self.data.loc[i, 'highway_str'] = 'track'
+                    else:
+                        self.data.loc[i, 'highway_str'] = 'mixed'
+                else:
+                    self.data.loc[i, 'highway_str'] = row["highway"]
+
+            # Unique list of road categories
+            self.road_categories = np.unique(self.data["highway_str"].values)
+            # Set every unknown road category to none
+            self.road_categories = [r if r in OsmGraph.road_colors else "none" for r in self.road_categories]
+
+        return(self)
 
     def plot(self, ax):
         """
         Plotting routine for OsmGraphs. Expects ax.
         """
-        if self.shape is None:
+        if self.data is None:
             print("! OsmGraph %s: Nothing to plot." % self.network)
         else:
             plot_kw = dict()
@@ -125,9 +163,29 @@ class OsmGraph:
                     color="k", alpha=0.3,
                     lw=1, ls=":")
 
-            self.shape.plot(
+            self.data.plot(
                 ax=ax, **plot_kw)
 
+    def plot_road_types(
+        self, ax,
+        legend_kw = dict(bbox_to_anchor=(1.01, 1), loc='upper left', frameon=False)
+        ):
+        """
+        Plotting routine for OsmGraphs road types. Expects ax.
+        """
+        if self.data is None or self.network != "roads":
+            print("! OsmGraph %s: Nothing to plot." % self.network)
+        else:
+            for cat in self.road_categories:
+                self.data.query("'%s' in highway_str" % cat).plot(
+                    ax=ax,
+                    label=cat,
+                    color=self.road_colors[cat],
+                    lw=1
+                    )
+    
+            ax.legend(**legend_kw)
+
 
 #######################################################
 ##### Automatic OSM Lookup from Overpass API

corny/figures.py

@@ -452,7 +452,7 @@ def add_circle(ax, x, y, radius=1,
 def add_colorbar(ax=None, points=None, label=None,
             ticks=None, ticklabels=None,
             ticks_kw=dict(),
-            bar_kw=dict(shrink=0.6, pad=0.02, aspect=20),
+            bar_kw=dict(shrink=0.6, pad=0.02, aspect=20, extend="both"),
             label_kw=dict(rotation=270, labelpad=20),
             ):
             
@@ -1413,7 +1413,8 @@ class Kriging:
                 self.z_max + self.z_resolution,
                 self.z_resolution),
                 vmin = self.z_min,
-                vmax = self.z_max
+                vmax = self.z_max,
+                extend="both"
             ))
         
     plot_z = plot_values
@@ -1472,7 +1473,22 @@ class Kriging:
         """
         Export to ESRI standard ASC grid file
         """
+        from corny.basics import latlon2utm
+        # Convert grids to UTM
+        y_grid, x = latlon2utm(lats=self.y_grid, lons=self.y_grid*0+self.x_grid[0])
+        y, x_grid = latlon2utm(lats=self.x_grid*0+self.y_grid[0], lons=self.x_grid)
+
+        # Fix regular grid spacing
+        x_step = np.mean(np.diff(x_grid))
+        for i in range(len(x_grid)):
+            x_grid[i] = x_grid[0] + i*x_step
+        y_step = np.mean(np.diff(y_grid))
+        for i in range(len(y_grid)):
+            y_grid[i] = y_grid[0] + i*y_step
+        
+        # Write ASCII file
+        from pykrige.kriging_tools import write_asc_grid
         write_asc_grid(
-            self.x_grid, self.y_grid,
+            x_grid, y_grid,
             self.z_grid,
-            filename=filename)
\ No newline at end of file
+            filename=filename)

default.cfg

@@ -601,7 +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
+# Special plots are: tubes, footprint, diurnal, map_xy, map_tiles, map_interpolation, roads
 PDF_plots = footprint
 # Displayed data range for vol. soil moisture in the plots.
 # Units: m³/m³

examples/rover.cfg

@@ -602,8 +602,8 @@ 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
+# Special plots are: tubes, footprint, diurnal, map_xy, map_tiles, map_interpolation, roads
+PDF_plots = footprint_length, footprint, tubes, map_grid, map_tiles, map_interpolation, roads
 # Displayed data range for vol. soil moisture in the plots.
 # Units: m³/m³
 PDF_sm_range = 0, 0.5
@@ -630,7 +630,7 @@ satellite_zoom = 16
 # - stamen           # Stamen terrain 
 # - toner            # Stamen toner
 # - watercolor       # Stamen water color
-tiles = satellite-google
+tiles = satellite-ms
 
 # Interpolation options if "map_interpolation" is in "PDF_plots".
 # Variogram model used for Ordinary Kriging interpolation (select one):

instantPASDy.py

@@ -1164,13 +1164,25 @@ def main(configfile=None, cfg_prefix=None):
             # back to Nfinal, but mavged
             data[Nfinal+'_temp_mov'] /= data['correct_bio']
             # N2SM first order
-            if config['conversion']['N0']:
-                N0 = config['conversion']['N0']
+            if cc(config,'conversion','N0'):
+                N0 = gc(config,'conversion','N0', dtype=float)
             else:
-                N0 = data[Nfinal].max()*1.1
+                N0 = np.nanmax(data[Nfinal])*1.1
                 print('! N0 is not defined, but we need it for road correction. Guessing %.0f cph.' % N0)
-            data['sm_temp'] = N2SM_Desilets(data[Nfinal+'_temp_mov'], N0, bd=data.bd, lw=data.lw, owe=data.owe,
-                        a0=float(config['conversion']['a0']), a1=float(config['conversion']['a1']), a2=float(config['conversion']['a2']))
+            
+            if cc(config,'conversion','N2sm_method','Desilets et al. (2010)'):    
+                data['sm_temp'] = N2SM_Desilets(data[Nfinal+'_temp_mov'], N0, bd=data.bd, lw=data.lw, owe=data.owe,
+                            a0=float(config['conversion']['a0']), a1=float(config['conversion']['a1']), a2=float(config['conversion']['a2']))
+            
+            elif cc(config,'conversion','N2sm_method','Schmidt et al. (2021)'):
+                if cc(config, 'conversion','Schmidt_parameterset'):
+                    parameterset = gc(config, 'conversion','Schmidt_parameterset')
+                else:
+                    parameterset = 'Mar21_mcnp_drf'
+                data['sm_temp'] = N2SM_Schmidt(data, Nfinal+'_temp_mov', ah, N0, bdstr='bd', lwstr='lw', owestr='owe', method=parameterset)
+            else:
+                print("! No valid sm conversion method selected.")    
+
             # C_road
             data['correct_road'] = 1/(1 + 0.42*(1-np.exp(-0.5*data.road_width)) * (1.06-4*data.road_moisture-(0.16 + data.road_moisture)/(0.39 + data['sm_temp'])))
             data.loc[(data['correct_road'] > 1), 'correct_road'] = 1
@@ -3664,6 +3676,64 @@ def main(configfile=None, cfg_prefix=None):
                             R.add_image(F.buff)
 
 
+                    if "roads" in pdf_plots_list:
+                        
+                        Pr.update("Roads")
+                        if not data[var].isnull().all() and lat != '' and lat in data.columns:
+
+                            from corny.figures import make_bbox
+                            from corny.exdata.openstreetmap import OsmGraph
+
+                            # Reduce data to the relevant pieces
+                            data_x = data[[lat,lon]].dropna()
+                            if "road" in data.columns:
+                                data_x["road"] = data["road"].fillna("none")
+
+                                # Set colors for road types
+                                data_x["road_color"] = [OsmGraph.road_colors[r] if r in OsmGraph.road_colors else "none" for r in data_x["road"]]
+                            else:
+                                data_x["road_color"] = "black"
+
+                            # Config
+                            resolution   = gc(config, "output", "interpolation_ll_resolution", alt=0.001, dtype=float)
+                            
+                            # Make bounding box from data
+                            bbox = make_bbox(
+                                data_x[lon].values,
+                                data_x[lat].values,
+                                pad=resolution)
+                            
+                            Roads  = OsmGraph("roads", bbox).get()
+
+                            R.add_page()
+                            R.add_title("Map of roads")
+                            R.add_par("The road network and road types have been downloaded from OpenStreetMap. Data within ~7 m received the road type attribute.")
+
+                            # Create the plot
+                            with Figure(
+                                size   = (9,9),
+                                # projection = "flat",
+                                extent = bbox,
+                                save="buff", #format="png"
+                                ) as F:
+
+                                ax = F.axes
+                                Roads.plot_road_types(ax=ax)
+
+                                data_x0 = data_x[data_x.road == "none"]
+                                data_x1 = data_x[data_x.road != "none"]
+                                ax.scatter(
+                                    data_x0[lon], data_x0[lat],
+                                    c=data_x0["road_color"], edgecolor="grey", lw=0.5, s=30,
+                                    zorder=2)
+                                ax.scatter(
+                                    data_x1[lon], data_x1[lat],
+                                    c=data_x1["road_color"], s=30,
+                                    zorder=3)
+                            
+                            R.add_image(F.buff)
+
+
                     for pdf_plots_var in pdf_plots_list:
 
                         if pdf_plots_var in ["footprint", "tubes", "diurnal"]: