In [1]:
import netCDF4
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
import pandas
import requests
import osgeo.osr
import tqdm
print(matplotlib.style.available)
matplotlib.style.use('bmh')
%matplotlib inline
In [20]:
# alternative (non-maintained? source https://data.overheid.nl/data/dataset/geotop)
url = 'http://dinodata.nl/opendap/GeoTOP/geotop.nc'
In [21]:
ds = netCDF4.Dataset(url)
In [73]:
x = ds.variables['x'][:]
y = ds.variables['y'][:]
z = ds.variables['z'][:]
z_0 = np.where(z == 0)[0]
z_min_10 = np.where(z == -10)[0]
litho = ds.variables['lithok'][::5, ::5, z_0]
In [74]:
_ = plt.hist(litho[~litho.mask])
In [75]:
labels = {
0: "mensen",
1: "veen",
2: "klei",
3: "zandige klei",
4: "lithoklasse 4",
5: "fijn zand",
6: "middelfijn zand",
7: "grof zand",
8: "grind",
9: "schelpen"
}
In [76]:
import matplotlib.colors
import matplotlib.cm
colors = [
'#3344ffff', # ophooglaag
'#845F4Cff', # organisch
'#734222ff', # klei
'#B99F71ff', # zandige klei
'#ff0000ff', # litho 4
'#E7D1C1ff', # fijn
'#c2b280ff', # middelfijn
'#969CAAff', # grof
'#D0D6D6ff', # grind
'#E5E5DBff' # schelp
]
geotop_cmap = matplotlib.colors.ListedColormap(
colors
)
In [77]:
fig, ax = plt.subplots(figsize=(13, 8))
ax.imshow(np.squeeze(litho.T), origin='top', extent=(x[0], x[-1], y[0], y[-1]), cmap=geotop_cmap)
for label in labels:
if label == 4:
continue
ax.plot(x[0], y[0], color=colors[label], label=labels[label])
ax.legend(loc='lower right')
Out[77]:
In [9]:
#
url = "https://waterwebservices.rijkswaterstaat.nl/METADATASERVICES_DBO/OphalenCatalogus/"
params = {
"CatalogusFilter":
{
"Eenheden": True,
"Grootheden":True,
"Hoedanigheden":True
}
}
resp = requests.post(url, json=params)
result = resp.json()
In [10]:
locaties = pandas.DataFrame.from_dict(result['LocatieLijst'])
locaties.head()
Out[10]:
In [11]:
aquo = pandas.DataFrame.from_dict(result['AquoMetadataLijst']).set_index('AquoMetadata_MessageID')
aquo.ix[60]
Out[11]:
In [12]:
aquolocaties = pandas.DataFrame.from_dict(result['AquoMetadataLocatieLijst'])
In [13]:
# only locations with waterlevel measurements relative to NAP
selected_locaties = aquolocaties[aquolocaties.AquoMetaData_MessageID == 60]
# get location information
selected_locaties = pandas.merge(locaties, selected_locaties)
# selected_locaties.X /= 10.0
# selected_locaties.Y /= 10.0
In [14]:
in_bbox = np.logical_and.reduce([
selected_locaties.X >= x[0],
selected_locaties.X < x[-1]
# selected_locaties.Y >= y[0],
# selected_locaties.Y < y[1],
])
selected_locaties.ix[in_bbox]
Out[14]:
In [15]:
plt.plot(selected_locaties.X, selected_locaties.Y, 'k.')
# plt.xlim(0.4e6, 0.9e6)
# plt.ylim(0.55e7, 0.60e7)
plt.xlim(0.4e5, 0.9e5)
plt.ylim(0.55e6, 0.60e6)
plt.plot(155000.0, 463000.0, 'rx')
Out[15]:
In [16]:
utm = osgeo.osr.SpatialReference()
rd = osgeo.osr.SpatialReference()
utm.ImportFromEPSG(25831)
# utm.ImportFromEPSG(32631)
rd.ImportFromEPSG(28992)
utm2rd = osgeo.osr.CoordinateTransformation(utm, rd)
pts = np.array(selected_locaties[['X', 'Y']])
pts_rd = np.array(utm2rd.TransformPoints(pts))
selected_locaties['X_rd'] = pts_rd[:, 0]
selected_locaties['Y_rd'] = pts_rd[:, 1]
In [17]:
fig, ax = plt.subplots(figsize=(13, 8))
ax.imshow(np.squeeze(litho.T), origin='top', extent=(x[0], x[-1], y[0], y[-1]), cmap=geotop_cmap)
ax.plot(selected_locaties.X_rd, selected_locaties.Y_rd, 'k.')
ax.set_xlim(x[0], x[-1])
ax.set_ylim(y[0], y[-1])
Out[17]:
In [18]:
selected_locaties[selected_locaties.Naam.apply(lambda x:'arl' in x)]
Out[18]:
In [19]:
stations = ['DENHDR', 'IJMH', 'HOEKVHLD', 'VLIS', 'VLISSGN', 'DELFZL', 'DLFZ', 'HARLGN', 'HARL']
In [20]:
main_stations = selected_locaties.ix[selected_locaties.Code.apply(lambda x: x in stations)].copy()
In [21]:
in_model = np.logical_and.reduce([
main_stations.X_rd >= x[0],
main_stations.X_rd < x[-1],
main_stations.Y_rd >= y[0],
main_stations.Y_rd < y[-1]
])
main_stations['in_model'] = in_model
main_stations['in_model'].all()
Out[21]:
In [22]:
main_stations['X_idx'] = np.searchsorted(x, main_stations.X_rd)
main_stations['Y_idx'] = np.searchsorted(y, main_stations.Y_rd)
In [23]:
lithos = {}
for i, station in main_stations.iterrows():
litho = ds.variables['lithok'][station.X_idx, station.Y_idx, :]
lithos[station.Code] = litho
In [24]:
litho_img = np.ma.masked_equal(list(lithos.values()), -127).T
fig, ax = plt.subplots(figsize=(13, 8))
ax.imshow(litho_img, origin='top', extent=(0, 9, z[0], z[-1]), cmap=geotop_cmap)
ax.set_aspect(0.05)
_ = ax.xaxis.set_ticks(list(range(1, len(stations) + 1)))
_ = ax.xaxis.set_ticklabels(stations, rotation=45)
In [25]:
fig, ax = plt.subplots()
ax.imshow(litho_img.copy())
ax.set_aspect(0.05)
In [26]:
litho_img.max()
Out[26]:
In [ ]: