diff --git a/examples/rover.cfg b/examples/rover.cfg
index 33fb798ab63c433df8bc6ab380b0ef9f96b282f9..202272b541aaebf73caf1b560573dd7dcaeff69c 100644
--- a/examples/rover.cfg
+++ b/examples/rover.cfg
@@ -95,6 +95,7 @@ missing_humidity = DWD
 missing_temperature = DWD
 # Range of corrected neutron counts in cph
 neutron_range = 3000, 20000
+neutron_range_per_tube = 1, 1000
 # Treat every change that is bigger than n_sigma * sqrt(N) as outlier
 neutron_n_sigma = 3
 # Range for other variables
@@ -251,7 +252,7 @@ measured_sm = 0.20
 ## Neutrons and soil moisture will be evaluated at that point
 ## Table syntax: name, lat, lon, radius (m)
 ## for multiline config make sure to add space before each line
-poi_table = 
+poi_table =
   river_plain,   51.877578, 12.244764,  50
   road_junction, 51.877191, 12.237821, 100
 ## How should distance weighting be performed?
@@ -271,7 +272,7 @@ make_CSV = yes
 ## which columns to save?
 ## e.g.: lat, lon, alt, pressure, temperature, relative_humidity, neutrons_raw, neutrons_proc, neutrons_proc_err, moisture_vol, moisture_vol_err_low, moisture_vol_err_upp, footprint_radius, footprint_depth
 ## Leave blank to save all columns
-CSV_columns = 
+CSV_columns =
 
 ## Date format in the output file (Y=year, m=month, d=day, H=hour, M=minute, S=second)
 CSV_datetime_format = Y-m-d H:M:S
@@ -287,7 +288,7 @@ CSV_NaN =
 # PDF
 make_PDF = yes
 ## Plot more columns (default is: N, pressure, rel humidity, temperature, abs humidity, moisture)
-PDF_plots = footprint_length, footprint, map_grid, map_tiles, map_grid_tiles, map_interpolation
+PDF_plots = footprint_length, footprint, tubes, map_grid, map_tiles_poi, map_grid, map_tiles, map_grid_tiles, map_interpolation
 ## Time format for x-axis in plots  (Y=year, m=month, d=day, H=hour, M=minute, S=second)
 PDF_time_format = H:M
 
@@ -314,5 +315,3 @@ KML_other_reverse = no
 ## Information on CRNS Footprint radius (in meter) and the representation of footprint geometry (point or polygons)
 make_footprint_polygons = yes
 # experimental: footprint_radius=50, footprint_type=polygon
-
-
diff --git a/instantPASDy.py b/instantPASDy.py
index 80b125fd0bac549ef427f83e64a76c3c0959521d..760c2618f512398f2e6826d3d76f0d27b44da273 100644
--- a/instantPASDy.py
+++ b/instantPASDy.py
@@ -1226,10 +1226,10 @@ def main(configfile=None):
                 poi_data.loc[i, smv+err]    = np.average(_data.loc[_query, smv+err],    weights=_data.loc[_query, 'distance_weight']) / np.sqrt(_data.loc[_query, smv].count())
                 poi_data.loc[i, smv+std]    =            _data.loc[_query, smv].std()
 
-                print('|   POI %15s (#%2.0f >%3.0f m): N = %5.0f +/-%4.0f (%4.0f) cph, SM = %4.1f +/-%4.1f (%4.1f) vol.%%' % (row['name'],
+                print('|   POI %15s (#%2.0f >%3.0f m): N = %5.0f +/-%4.0f cph, SM = %4.1f +/-%4.1f vol.%%' % (row['name'],
                     _data.loc[_query, 'distance2poi'].count(), _data.loc[_query, 'distance2poi'].min(),
-                    poi_data.loc[i, Nfinal],  poi_data.loc[i, Nfinal+err],  poi_data.loc[i, Nfinal+std],
-                    poi_data.loc[i, smv]*100, poi_data.loc[i, smv+err]*100, poi_data.loc[i, smv+std]*100))
+                    poi_data.loc[i, Nfinal],  poi_data.loc[i, Nfinal+err], # poi_data.loc[i, Nfinal+std],
+                    poi_data.loc[i, smv]*100, poi_data.loc[i, smv+err]*100)) #, poi_data.loc[i, smv+std]*100)) # (%4.1f)
             else:
                 print('|   POI %15s (  >%5.0f m): (not covered)' % (row['name'], _data.distance2poi.min()))
 
@@ -1586,9 +1586,10 @@ def main(configfile=None):
                                 ax.grid(b=True, color='black', alpha=0.2)
                                 ax.text(data[var].index[0], ymax*0.95, var)
 
-                        for i in range(int(mysize)*int(mysize)-len(N_tubes)):
-                            with Fig(ylabel='N', ylim=(0,3)) as ax:
-                                ax.set_axis_off()
+                        if mylayout != (2,1) and len(N_tubes)>=2:
+                            for i in range(int(mysize)*int(mysize)-len(N_tubes)):
+                                with Fig(ylabel='N', ylim=(0,3)) as ax:
+                                    ax.set_axis_off()
                             #with Fig(ylabel='N', ylim=(0,3)):
                             #    pass #plt.plot(data[N], color='black', ls='', marker='o', ms=2, markeredgewidth=0)
 
@@ -1728,7 +1729,104 @@ def main(configfile=None):
 
                         if cc(config,'conversion','poi_plot','yes'):
                             for i, poi in poi_data.iterrows():
-                                ax.scatter( poi.lon, poi.lat, c='black', marker='+', s=100)
+                                ax.scatter( poi.lon, poi.lat, c='black', marker='+', s=100, transform=ccrs.PlateCarree())
+
+
+                ## Plot with tiles (EXPERIMENTAL)
+                elif var == 'map_tiles_poi':
+
+                    var = sm
+                    # TODO Those params should go to the config file
+                    todo_config_zoom = 17
+                    todo_config_map = 'satellite'
+                    todo_config_markersize = 15
+
+                    print('POI plot with tiles... ', end='')
+                    datall = data.dropna(subset=[sm]).copy()
+
+                    if todo_config_map == 'satellite':
+                        cimgt.QuadtreeTiles.get_image = image_spoof # reformat web request for street map spoofing
+                        osm_img = cimgt.QuadtreeTiles() # spoofed, downloaded street map
+                    elif todo_config_map == 'OSM':
+                        cimgt.OSM.get_image = image_spoof # reformat web request for street map spoofing
+                        osm_img = cimgt.OSM() # spoofed, downloaded street map
+                    else:
+                        todo_config_map = None
+
+                    lonpad = 0.1*(datall[lon].max()-datall[lon].min())
+                    latpad = 0.3*(datall[lat].max()-datall[lat].min())
+
+                    datax = data[[lon,lat,var]].dropna()
+                    if len(datax)==0:
+                        raise ValueError('No records with valid sm left after processing, stop printing')
+
+                    ## Deprecated:
+                    # import pyproj
+                    # xs, ys = pyproj.transform(pyproj.Proj("epsg:4326"), pyproj.Proj("epsg:31468"), datax[lat].values, datax[lon].values)
+                    ## New pyproj2 syntax:
+                    from pyproj import Transformer
+                    xs, ys = Transformer.from_crs("epsg:4326", "epsg:31468").transform(datax[lat].values, datax[lon].values)
+
+                    datax['y'] = xs
+                    datax['x'] = ys
+
+                    from scipy.interpolate import griddata
+                    # data coordinates and values
+                    x = datax.x.values - datax.x.min()
+                    y = datax.y.values - datax.y.min()
+                    z = datax[var].values
+
+                    smrange = (0,0.5)
+
+                    xrange = x.max()-x.min()
+                    yrange = y.max()-y.min()
+                    if xrange==0: xrange=1
+                    if yrange==0: yrange=1
+                    xpad = 0.05*xrange
+                    ypad = 0.10*yrange
+
+                    mysize = (10,10/xrange*yrange)
+                    if mysize[1] > 20:
+                        mysize = (10*xrange/yrange,10)
+
+                    zoom = todo_config_zoom #18-int(15*(data[lon].max()-data[lon].min())) # 0->18, 0.5->16, 1->14
+                    # Do not choose >=17 for more than a small field site! Python will crash.
+                    #print("(zoom %i) " % zoom, end='')
+
+                    mygrids = np.ceil(np.sqrt(len(poi_data)))
+                    if len(poi_data) == 2:
+                        mylayout = (2,1)
+                    else:
+                        mylayout = (int(mygrids),int(mygrids))
+
+                    with Fig(fig, layout=mylayout) as axes:
+                        for i, poi in poi_data.iterrows():
+                            xlimits = (poi.lon-0.0014, poi.lon+0.0014)
+                            ylimits = (poi.lat-0.0008, poi.lat+0.0008)
+
+                            with Fig(title='%s (100 m radius)' % poi['name'], xlabel='lon', ylabel='lat', time_series=False, proj=ccrs.PlateCarree()) as ax:
+                                ax.scatter( datall[lon], datall[lat], c='black', transform=ccrs.PlateCarree(), s=todo_config_markersize*2)
+                                q = ax.scatter( datall[lon], datall[lat], c=datall[sm], transform=ccrs.PlateCarree(),
+                                                s=todo_config_markersize, cmap='Spectral', vmin=smrange[0], vmax=smrange[1])
+                                ax.set_xlim(xlimits)
+                                ax.set_ylim(ylimits)
+                                #ax.set_extent([xlimits[1], ylimits[1], xlimits[0], ylimits[0]], crs=ccrs.PlateCarree())
+                                if not todo_config_map is None:
+                                    ax.add_image(osm_img, zoom)
+
+                                #cb = plt.colorbar(q, ax=ax, shrink=0.7, pad=0.03, aspect=30)
+                                #cb.set_label("vol. soil moisture (in m$^3$/m$^3$)", rotation=270, labelpad=20)
+
+                                #if cc(config,'conversion','poi_plot','yes'):
+                                #    for i, poi in poi_data.iterrows():
+                                ax.scatter( poi.lon, poi.lat, c='black', marker='+', s=100, transform=ccrs.PlateCarree())
+
+                        if mylayout != (2,1):
+                            for i in range(int(mygrids)*int(mygrids)-len(poi_data)):
+                                with Fig(ylabel='N', ylim=(0,3)) as ax:
+                                    ax.set_axis_off()
+
+
 
                 ## Plot grid with tiles (EXPERIMENTAL)
                 elif var == 'map_grid_tiles':
@@ -1875,6 +1973,7 @@ def main(configfile=None):
                             for i, poi in poi_data.iterrows():
                                 ax.scatter( poi.utmx-datax.x.min(), poi.utmy-datax.y.min(), c='black', marker='+', s=100)
 
+
                 elif var in data.columns:
                     with Fig(fig, title=var, ylabel=var):
                         if not data[var].isnull().all():