Commit 35a65276 authored by Martin Schrön's avatar Martin Schrön
Browse files

added POI map plot feature and tube plot for example rover config

parent 25dca0ef
......@@ -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
......@@ -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():
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment