Merge pull request #361 from gwmod/improve_geotop_strat_units

Add information about geotop stratigraphic units

nlmod/data/geotop/REF_GTP_LITHO_CLASS.csv

@@ -0,0 +1,12 @@
+LITHO_CLASS_CD,DESCRIPTION,VOXEL_NR,SEQ_NR,RED_DEC,GREEN_DEC,BLUE_DEC
+a,antropogeen,0,0,200,200,200
+v,organisch materiaal (veen),1,1,157,78,64
+k,klei,2,2,0,146,0
+kz,"kleiig zand, zandige klei en leem",3,3,194,207,92
+zf,zand fijn,5,5,255,255,0
+zm,zand midden,6,6,243,225,6
+zg,zand grof,7,7,231,195,22
+g,grind,8,8,216,163,32
+she,schelpen,9,9,95,95,255
+z,zand niet gedifferentieerd,Null,Null,243,225,6
+o,overig,Null,Null,144,144,144

nlmod/data/geotop/REF_GTP_STR_UNIT.csv

@@ -0,0 +1,107 @@
+STR_UNIT_CD,DESCRIPTION,VOXEL_NR,SEQ_NR,RED_DEC,GREEN_DEC,BLUE_DEC
+AAOP,"Antropogene afzettingen, opgebrachte grond",1000,1,200,200,200
+AAES,"Antropogene afzettingen, esdekken",1005,2,110,110,110
+NIGR,"Formatie van Nieuwkoop, Laagpakket van Griendtsveen",1010,3,132,88,44
+NINB,"Formatie van Nieuwkoop, Laag van Nij Beets",1045,4,180,90,20
+NASC,"Formatie van Naaldwijk, Laagpakket van Schoorl",1020,5,255,255,160
+ONAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren (gedeelte boven NAZA)",1030,6,152,220,117
+NAZA,"Formatie van Naaldwijk, Laagpakket van Zandvoort",1040,7,254,183,56
+NAWAZU,"Formatie van Naaldwijk, Laagpakket van Walcheren, Zuiderzee Laag",1048,14,177,160,199
+NAWAAL,"Formatie van Naaldwijk, Laagpakket van Walcheren, Almere Laag",1049,15,96,73,122
+NAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren, gelegen onder Formatie van Naaldwijk, Laagpakket van Zandvoort",1050,16,152,220,117
+BHEC,Formatie van Echteld (gedeelte buiten NIHO),1060,,170,255,245
+OEC,Formatie van Echteld (gedeelte boven NIHO),1070,13,170,255,245
+NAWOBE,"Formatie van Naaldwijk, Laagpakket van Wormer, Laag van Bergen",1080,18,0,157,155
+KK1,Kreekrak Formatie (gedeelte boven NIHO),1085,17,68,138,112
+NIFL,"Formatie van Nieuwkoop, Flevomeer Laag",1089,19,252,213,180
+NIHO,"Formatie van Nieuwkoop, Hollandveen Laagpakket",1090,20,208,130,40
+NAZA2,"Formatie van Naaldwijk, Laagpakket van Zandvoort (gedeelte onder NIHO)",1095,21,254,183,56
+NAWO,"Formatie van Naaldwijk, Laagpakket van Wormer",1100,27,63,165,76
+NWNZ,"Formatie van Naaldwijk, laagpakketten van Wormer en Zandvoort",1110,,63,165,76
+NAWOVE,"Formatie van Naaldwijk, Laagpakket van Wormer, Laag van Velsen",1120,29,189,183,107
+KK2,Kreekrak Formatie (gedeelte onder NIHO),1125,30,68,138,112
+NIBA,"Formatie van Nieuwkoop, Basisveen Laag",1130,31,152,47,10
+NA,Formatie van Naaldwijk,2000,12,105,200,100
+EC,Formatie van Echteld,2010,28,170,255,245
+NI,Formatie van Nieuwkoop,2020,,208,130,40
+KK,Kreekrak Formatie,2030,,68,138,112
+BXKO,"Formatie van Boxtel, Laagpakket van Kootwijk",3000,34,255,255,190
+BXSI,"Formatie van Boxtel, Laagpakket van Singraven",3010,,190,190,0
+BXSI1,Formatie van Boxtel Laagpakket van Singraven (bovenste deel),3011,35,190,190,0
+BXWI,"Formatie van Boxtel, Laagpakket van Wierden",3020,36,255,255,80
+BXSI2,Formatie van Boxtel Laagpakket van Singraven (onderste deel),3012,37,255,255,130
+BXWIKO,"Formatie van Boxtel, laagpakketten van Wierden en Kootwijk",3025,38,255,255,80
+BXWISIKO,"Formatie van Boxtel, laagpakketten van Wierden, Singraven en Kootwijk",3030,39,255,255,80
+BXDE,"Formatie van Boxtel, Laagpakket van Delwijnen",3040,40,225,225,130
+BXDEKO,"Formatie van Boxtel, laagpakketten van Delwijnen en Kootwijk",3045,41,225,225,130
+BXSC,"Formatie van Boxtel, Laagpakket van Schimmert",3050,42,255,205,0
+BXLM,"Formatie van Boxtel, Laagpakket van Liempde",3060,43,225,190,0
+BXBS,"Formatie van Boxtel, Laagpakket van Best",3090,46,235,205,0
+BX,Formatie van Boxtel,3100,47,255,235,0
+KRWY,"Formatie van Kreftenheye, Laag van Wijchen",4000,48,86,0,0
+KRBXDE,"Formatie van Kreftenheye en Formatie van Boxtel, Laagpakket van Delwijnen",4010,49,176,48,96
+KRZU,"Formatie van Kreftenheye, Laagpakket van Zutphen",4020,50,136,0,0
+KROE,"Formatie van Kreftenheye, gelegen onder de Eem Formatie",4030,61,176,48,96
+KRTW,"Formatie van Kreftenheye, Laagpakket van Twello",4040,51,156,18,46
+KR,Formatie van Kreftenheye,4050,52,176,48,96
+BEOM,"Formatie van Beegden, Laagpakket van Oost-Maarland",4055,32,129,105,123
+BEWY,"Formatie van Beegden, Laag van Wijchen",4060,53,170,155,180
+BERO,"Formatie van Beegden, Laag van Rosmalen",4070,54,160,140,155
+BE,Formatie van Beegden,4080,55,200,200,255
+KW1,Formatie van Koewacht (kleiige top),4085,56,152,139,0
+KW,Formatie van Koewacht,4090,57,172,169,43
+WB,Formatie van Woudenberg,4100,58,137,67,30
+EE,Eem Formatie,4110,59,210,255,115
+EEWB,Formatie van Woudenberg en Eem Formatie,4120,60,210,255,115
+DR,Formatie van Drente,5000,62,255,127,80
+DRGI,"Formatie van Drente, Laagpakket van Gieten",5010,63,235,97,30
+GE,Door landijs gestuwde afzettingen,5020,64,156,156,156
+DN,Formatie van Drachten,5030,65,250,250,210
+URTY,"Formatie van Urk, Laagpakket van Tynje",5040,66,169,163,87
+PE,Formatie van Peelo,5050,67,238,130,238
+UR,Formatie van Urk,5060,68,189,183,107
+ST,Formatie van Sterksel,5070,69,205,92,92
+AP,Formatie van Appelscha,5080,70,218,165,32
+SY,Formatie van Stramproy,5090,71,255,228,181
+PZ,Formatie van Peize,5100,72,255,255,0
+WA,Formatie van Waalre,5110,73,255,165,0
+PZWA,Formatie van Peize en Formatie van Waalre,5120,74,255,204,0
+MS,Formatie van Maassluis,5130,75,135,206,235
+KI,Kiezeloöliet Formatie,5140,76,188,143,143
+OO,Formatie van Oosterhout,5150,77,118,157,39
+IE,Formatie van Inden,5160,78,236,121,193
+VI,Formatie van Ville,5170,79,153,102,0
+BR,Formatie van Breda,5180,80,108,188,150
+VE,Formatie van Veldhoven,5185,,102,100,16
+RUBO,"Rupel Formatie, Laagpakket van Boom",5190,81,154,78,163
+RU,Rupel Formatie,5200,82,184,123,238
+TOZEWA,"Formatie van Tongeren, Laagpakket van Zelzate, Laag van Watervliet",5210,83,80,144,194
+TOGO,"Formatie van Tongeren, Laagpakket van Goudsberg",5220,84,70,129,169
+TO,Formatie van Tongeren,5230,85,90,159,219
+DOAS,"Formatie van Dongen, Laagpakket van Asse",5240,86,186,146,141
+DOIE,"Formatie van Dongen, Laagpakket van Ieper",5250,87,206,176,191
+DO,Formatie van Dongen,5260,88,216,191,216
+LA,Formatie van Landen,5270,89,208,32,144
+HT,Formatie van Heijenrath,5280,90,178,34,34
+HO,Formatie van Houthem,5290,91,210,105,30
+MT,Formatie van Maastricht,5300,92,255,160,102
+GU,Formatie van Gulpen,5310,93,245,222,179
+VA,Formatie van Vaals,5320,94,21,153,79
+AK,Formatie van Aken,5330,95,152,231,205
+AEC,Formatie van Echteld (geulafzettingen generatie A),6000,9,118,147,60
+ABEOM,"Formatie van Beegden, Laagpakket van Oost-Maarland (geulafzettingen generatie A)",6005,33,118,147,60
+ANAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren (geulafzettingen generatie A)",6010,8,118,147,60
+ANAWO,"Formatie van Naaldwijk, Laagpakket van Wormer (geulafzettingen generatie A)",6020,,118,147,60
+BEC,Formatie van Echteld (geulafzettingen generatie B),6100,11,102,205,171
+BNAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren (geulafzettingen generatie B)",6110,10,102,205,171
+BNAWO,"Formatie van Naaldwijk, Laagpakket van Wormer (geulafzettingen generatie B)",6120,,102,205,171
+CEC,Formatie van Echteld (geulafzettingen generatie C),6200,22,170,196,255
+CNAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren (geulafzettingen generatie C)",6210,,170,196,255
+CNAWO,"Formatie van Naaldwijk, Laagpakket van Wormer (geulafzettingen generatie C)",6220,,170,196,255
+DEC,Formatie van Echteld (geulafzettingen generatie D),6300,24,102,153,205
+DNAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren (geulafzettingen generatie D)",6310,,102,153,205
+DNAWO,"Formatie van Naaldwijk, Laagpakket van Wormer (geulafzettingen generatie D)",6320,23,102,153,205
+EEC,Formatie van Echteld (geulafzettingen generatie E),6400,26,27,101,175
+ENAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren (geulafzettingen generatie E)",6410,,27,101,175
+ENAWO,"Formatie van Naaldwijk, Laagpakket van Wormer (geulafzettingen generatie E)",6420,25,27,101,175
+NN,Niet formeel ingedeelde afzettingen of onbekend,0,,85,85,85

nlmod/data/geotop/geo_eenheden.csv

@@ -6,13 +6,14 @@ strat,code,name
 1030,ONAWA,"Formatie van Naaldwijk, Laagpakket van
 Walcheren, gelegen boven Formatie van Naaldwijk, Laagpakket van Zandvoort"
 1040,NAZA,"Formatie van Naaldwijk, Laagpakket van Zandvoort"
+1045,NINB,"Formatie van Nieuwkoop, Laag van Nij Beets"
 1050,NAWA,"Formatie van Naaldwijk, Laagpakket van Walcheren, gelegen onder Formatie van Naaldwijk, Laagpakket van Zandvoort"
 1060,BHEC,"Formatie van Echteld, gelegen buiten de
 verbreiding van de Formatie van Nieuwkoop, Hollandveen Laagpakket"
 1070,OEC,"Formatie van Echteld, gelegen boven de Formatie van Nieuwkoop, Hollandveen Laagpakket"
 1080,NAWOBE,"Formatie van Naaldwijk, Laagpakket van Wormer, Laag van Bergen"
 1090,NIHO,"Formatie van Nieuwkoop, Hollandveen Laagpakket"
-1095,NINB,"Formatie van Nieuwkoop, Laagpakket van Nij Beets"
+1095,NAZA2,"Formatie van Naaldwijk, Laagpakket van Zandvoort (gedeelte onder NIHO)"
 1100,NAWO,"Formatie van Naaldwijk, Laagpakket van Wormer"
 1110,NWNZ,"Formatie van Naaldwijk, Laagpakketten van
 Wormer en Zandvoort (gecombineerde eenheid in modelgebied Zeeland)"

nlmod/dims/base.py

@@ -206,7 +206,7 @@ def to_model_ds(
  return ds
 
 
-def extrapolate_ds(ds, mask=None, layer="layer"):
+def extrapolate_ds(ds, mask=None, layer="layer", mask_values=None):
  """Fill missing data in layermodel based on nearest interpolation.
 
  Used for ensuring layer model contains data everywhere. Useful for
@@ -223,6 +223,10 @@ def extrapolate_ds(ds, mask=None, layer="layer"):
  variable. The default is None.
  layer : str, optional
  The name of the layer dimension. The default is 'layer'.
+ mask_values : np.ndarray, optional
+ Boolean mask for each cell, with a value of True if its value is used to fill
+ data in mask. When mask_values is None, it is determined from mask. The default
+ is None.
 
  Returns
  -------
@@ -234,6 +238,8 @@ def extrapolate_ds(ds, mask=None, layer="layer"):
  if not mask.any():
  # all of the model cells are is inside the known area
  return ds
+ if mask_values is None:
+ mask_values = ~mask
  if mask.all():
  raise (ValueError("The model only contains NaNs"))
  if "gridtype" in ds.attrs and ds.gridtype == "vertex":
@@ -243,7 +249,7 @@ def extrapolate_ds(ds, mask=None, layer="layer"):
  else:
  x, y = np.meshgrid(ds.x, ds.y)
  dims = ("y", "x")
- points = np.stack((x[~mask], y[~mask]), axis=1)
+ points = np.stack((x[mask_values], y[mask_values]), axis=1)
  xi = np.stack((x[mask], y[mask]), axis=1)
  # geneterate the tree only once, to increase speed
  tree = cKDTree(points)
@@ -254,11 +260,11 @@ def extrapolate_ds(ds, mask=None, layer="layer"):
  data = ds[key].data
  if ds[key].dims == dims:
  if np.isnan(data[mask]).sum() > 0: # do not update if no NaNs
- data[mask] = data[~mask][i]
+ data[mask] = data[mask_values][i]
  elif ds[key].dims == (layer,) + dims:
  for lay in range(len(ds[layer])):
  if np.isnan(data[lay, mask]).sum() > 0: # do not update if no NaNs
- data[lay, mask] = data[lay, ~mask][i]
+ data[lay, mask] = data[lay, mask_values][i]
  else:
  logger.warning(
  f"Data variable '{key}' not extrapolated because "

nlmod/dims/grid.py

@@ -702,8 +702,6 @@ def update_ds_from_layer_ds(ds, layer_ds, method="nearest", **kwargs):
  Dataset with variables from layer_ds.
  """
  if not layer_ds.layer.equals(ds.layer):
- # do not change the original Dataset
- layer_ds = layer_ds.copy()
  # update layers in ds
  drop_vars = []
  for var in ds.data_vars:
@@ -1201,7 +1199,7 @@ def gdf_to_data_array_struc(
  DESCRIPTION.
  """
  warnings.warn(
- "The method gdf_to_data_array_struc is deprecated. Please use gdf_to_da instead",
+ "The method gdf_to_data_array_struc is deprecated. Please use gdf_to_da instead.",
  DeprecationWarning,
  )
 
@@ -1770,7 +1768,7 @@ def get_thickness_from_topbot(top, bot):
  or (layer, icell2d).
  """
  warnings.warn(
- "The method get_thickness_from_topbot is deprecated. Please use nlmod.layers.calculate_thickness instead",
+ "The method get_thickness_from_topbot is deprecated. Please use nlmod.layers.calculate_thickness instead.",
  DeprecationWarning,
  )
 

nlmod/dims/layers.py

@@ -472,6 +472,7 @@ def kheq_combined_layers(kh, thickness, reindexer):
  kheq = np.nansum(
  thickness.data[v, :, :] * kh.data[v, :, :], axis=0
  ) / np.nansum(thickness.data[v, :, :], axis=0)
+ kheq[np.isinf(kheq)] = np.nan
  else:
  kheq = kh.data[v]
  da_kh.data[k] = kheq
@@ -512,6 +513,7 @@ def kveq_combined_layers(kv, thickness, reindexer):
  kveq = np.nansum(thickness.data[v, :, :], axis=0) / np.nansum(
  thickness.data[v, :, :] / kv.data[v, :, :], axis=0
  )
+ kveq[np.isinf(kveq)] = np.nan
  else:
  kveq = kv.data[v]
  da_kv.data[k] = kveq
@@ -536,7 +538,7 @@ def combine_layers_ds(
  ds : xarray.Dataset
  xarray Dataset containing information about layers
  (layers, top and bot)
- combine_layers : list of tuple of ints
+ combine_layers : list of tuple of ints, or dict of layer names
  list of tuples, with each tuple containing integers indicating
  layer indices to combine into one layer. E.g. [(0, 1), (2, 3)] will
  combine layers 0 and 1 into a single layer (with index 0) and layers
@@ -586,6 +588,17 @@ def combine_layers_ds(
  ds = ds.copy()
  ds["top"] = ds["botm"] + calculate_thickness(ds)
 
+ if isinstance(combine_layers, dict):
+ new_layer_names = combine_layers.keys()
+ combine_layers = [
+ tuple(np.where(ds.layer.isin(x))[0]) for x in combine_layers.values()
+ ]
+ # make sure there are no layers in between
+ assert np.all([(np.diff(x) == 1).all() for x in combine_layers])
+ else:
+ new_layer_names = [ds.layer.data[x[0]] for x in combine_layers]
+ new_layer_names = dict(zip(combine_layers, new_layer_names))
+
  da_dict = {}
 
  new_top, new_bot, reindexer = layer_combine_top_bot(
@@ -604,19 +617,20 @@ def combine_layers_ds(
  if kv is not None:
  logger.info(f"Calculate equivalent '{kv}' for combined layers.")
  da_dict[kv] = kveq_combined_layers(ds[kv], thickness, reindexer)
- if kD is not None:
+ if kD is not None and kD in ds:
  logger.info(f"Calculate value '{kD}' for combined layers with sum.")
  da_dict[kD] = sum_param_combined_layers(ds[kD], reindexer)
- if c is not None:
+ if c is not None and c in ds:
  logger.info(f"Calculate value '{c}' for combined layers with sum.")
  da_dict[c] = sum_param_combined_layers(ds[c], reindexer)
 
  # get new layer names, based on first sub-layer from each combined layer
  layer_names = []
  for _, i in reindexer.items():
  if isinstance(i, tuple):
- i = i[0]
- layercode = ds[layer].data[i]
+ layercode = new_layer_names[i]
+ else:
+ layercode = ds[layer].data[i]
  layer_names.append(layercode)
 
  # assign new layer names

nlmod/dims/resample.py

@@ -583,7 +583,7 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.nan):
 def extent_to_polygon(extent):
  logger.warning(
  "nlmod.resample.extent_to_polygon is deprecated. "
- "Use nlmod.util.extent_to_polygon instead"
+ "Use nlmod.util.extent_to_polygon instead."
  )
  from ..util import extent_to_polygon
 
@@ -594,7 +594,7 @@ def get_extent_polygon(ds, rotated=True):
  """Get the model extent, as a shapely Polygon."""
  logger.warning(
  "nlmod.resample.get_extent_polygon is deprecated. "
- "Use nlmod.grid.get_extent_polygon instead"
+ "Use nlmod.grid.get_extent_polygon instead."
  )
  from .grid import get_extent_polygon
 
@@ -605,7 +605,7 @@ def affine_transform_gdf(gdf, affine):
  """Apply an affine transformation to a geopandas GeoDataFrame."""
  logger.warning(
  "nlmod.resample.affine_transform_gdf is deprecated. "
- "Use nlmod.grid.affine_transform_gdf instead"
+ "Use nlmod.grid.affine_transform_gdf instead."
  )
  from .grid import affine_transform_gdf
 
@@ -615,7 +615,7 @@ def affine_transform_gdf(gdf, affine):
 def get_extent(ds, rotated=True):
  """Get the model extent, corrected for angrot if necessary."""
  logger.warning(
- "nlmod.resample.get_extent is deprecated. Use nlmod.grid.get_extent instead"
+ "nlmod.resample.get_extent is deprecated. Use nlmod.grid.get_extent instead."
  )
  from .grid import get_extent
 
@@ -626,7 +626,7 @@ def get_affine_mod_to_world(ds):
  """Get the affine-transformation from model to real-world coordinates."""
  logger.warning(
  "nlmod.resample.get_affine_mod_to_world is deprecated. "
- "Use nlmod.grid.get_affine_mod_to_world instead"
+ "Use nlmod.grid.get_affine_mod_to_world instead."
  )
  from .grid import get_affine_mod_to_world
 
@@ -637,7 +637,7 @@ def get_affine_world_to_mod(ds):
  """Get the affine-transformation from real-world to model coordinates."""
  logger.warning(
  "nlmod.resample.get_affine_world_to_mod is deprecated. "
- "Use nlmod.grid.get_affine_world_to_mod instead"
+ "Use nlmod.grid.get_affine_world_to_mod instead."
  )
  from .grid import get_affine_world_to_mod
 
@@ -647,7 +647,7 @@ def get_affine_world_to_mod(ds):
 def get_affine(ds, sx=None, sy=None):
  """Get the affine-transformation, from pixel to real-world coordinates."""
  logger.warning(
- "nlmod.resample.get_affine is deprecated. Use nlmod.grid.get_affine instead"
+ "nlmod.resample.get_affine is deprecated. Use nlmod.grid.get_affine instead."
  )
  from .grid import get_affine