213 fixes related to xugrid issues

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.11.6
+current_version = 0.11.7
 commit = True
 tag = True
 

CITATION.cff

@@ -8,5 +8,5 @@ authors:
     orcid: https://orcid.org/0000-0002-6349-818X
 title: "dfm_tools: A Python package for pre- and postprocessing D-FlowFM model input and output files"
 type: software
-version: 0.11.6
+version: 0.11.7
 doi: 

dfm_tools/__init__.py

@@ -4,7 +4,7 @@
 
 __author__ = """Jelmer Veenstra"""
 __email__ = '[email protected]'
-__version__ = '0.11.6'
+__version__ = '0.11.7'
 
 from dfm_tools.deprecated import *
 from dfm_tools.errors import *

dfm_tools/get_nc.py

@@ -472,37 +472,30 @@ def rasterize_ugrid(uds:xu.UgridDataset, ds_like:xr.Dataset = None, resolution:f
         DESCRIPTION.
 
     """
-    #TODO: maybe put part of code in xugrid (https://github.com/Deltares/xugrid/issues/31)
-    #TODO: vars can also be rasterized with uds_facevars[var].ugrid.rasterize(resolution), but is not efficient. Wait for uds.rasterize() method: https://github.com/Deltares/xugrid/issues/61
-    if not isinstance(uds,xu.core.wrap.UgridDataset):
-        raise TypeError(f'rasterize_ugrid expected xu.core.wrap.UgridDataset, got {type(uds)} instead')
 
-    grid = uds.grid
-    xu_facedim = uds.grid.face_dimension
-    uds_facevars = Dataset_varswithdim(uds,xu_facedim)
-
-    if ds_like is not None:
-        regx = ds_like.x
-        regy = ds_like.y
+    if isinstance(uds,xu.core.wrap.UgridDataset):
+        xu_facedim = uds.grid.face_dimension
+        uds_facevars = Dataset_varswithdim(uds,xu_facedim)
+        face_str = f'Dataset with {len(uds_facevars.data_vars)} face variables'
+    elif isinstance(uds,xu.core.wrap.UgridDataArray):
+        face_str = 'DataArray'
     else:
-        xmin, ymin, xmax, ymax = grid.bounds
+        raise TypeError(f'rasterize_ugrid expected xu.core.wrap.UgridDataset or xu.core.wrap.UgridDataArray, got {type(uds)} instead')
+
+    if ds_like is None:
+        xmin, ymin, xmax, ymax = uds.grid.bounds
         dx = xmax - xmin
         dy = ymax - ymin
         if resolution is None: # check if a rasterization resolution is passed, otherwise default to 200 raster cells otherwise for the smallest axis.
             resolution = min(dx, dy) / 200
         d = abs(resolution)
         regx = np.arange(xmin + 0.5 * d, xmax, d)
         regy = np.arange(ymin + 0.5 * d, ymax, d)
+        ds_like = xr.DataArray(np.empty((len(regy), len(regx))), {"y": regy, "x": regx}, ["y", "x"])
 
-    regx, regy, index = grid.rasterize_like(x=regx,y=regy) #TODO: this can be used to steer rasterization, eg with xstart/ystart/xres/yres
-    index_da = xr.DataArray(index,dims=('y','x'))
-
-    print(f'>> rasterizing ugrid dataset with {len(uds_facevars.data_vars)} face variables to shape={index_da.shape}: ',end='')
+    print(f'>> rasterizing ugrid {face_str} to shape=({len(ds_like.y)},{len(ds_like.y)}): ',end='')
     dtstart = dt.datetime.now()
-    ds = uds_facevars.isel({xu_facedim:index_da})
-    ds = ds.where(index_da != grid.fill_value)
-    ds['x'] = xr.DataArray(regx,dims='x')
-    ds['y'] = xr.DataArray(regy,dims='y')
+    ds = uds.ugrid.rasterize_like(other=ds_like)
     print(f'{(dt.datetime.now()-dtstart).total_seconds():.2f} sec')
 
     return ds

environment.yml

@@ -20,7 +20,7 @@ dependencies:
   - dask
   - netcdf4>=1.5.3
   - bottleneck
-  - xugrid>=0.2.1
+  - xugrid>=0.4.0
   - cdsapi
   - pydap>=3.3.0
   - pip:

requirements.txt

@@ -11,7 +11,7 @@ xarray<2023.3.0
 dask
 netcdf4>=1.5.3
 bottleneck
-xugrid>=0.2.1
+xugrid>=0.4.0
 cdsapi
 pydap>=3.3.0
 hydrolib-core>=0.5.1

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = dfm_tools
-version = 0.11.6
+version = 0.11.7
 author = Jelmer Veenstra
 author_email = [email protected]
 description = dfm_tools are pre- and post-processing tools for Delft3D FM
@@ -32,7 +32,7 @@ install_requires =
 	dask
 	netcdf4>=1.5.3
 	bottleneck
-	xugrid>=0.2.1
+	xugrid>=0.4.0
 	cdsapi
 	pydap>=3.3.0
 	hydrolib-core>=0.5.1

tests/examples/postprocess_mapnc_regular.py

@@ -78,7 +78,7 @@
         data_u_thin = data_u_tsel.loc[::thinning,::thinning]
         data_v_thin = data_v_tsel.loc[::thinning,::thinning]
         ax.quiver(data_u_thin[name_uv_x], data_u_thin[name_uv_y], data_u_thin, data_v_thin, #raises warnings for some reason
-                  color='w')#,scale=50,width=0.008)#, edgecolor='face', cmap='jet')
+                  color='w')#,scale=50,width=0.008), cmap='jet')
     fig.tight_layout()
     plt.savefig(os.path.join(dir_output,f'{basename}_magn_pcolorquiver'))
 

tests/examples/postprocess_mapnc_ugrid.py

@@ -20,7 +20,7 @@
 file_nc_list = [os.path.join(dir_testinput,'DFM_curvedbend_3D','cb_3d_map.nc'), #sigmalayer
                 os.path.join(dir_testinput,'DFM_grevelingen_3D','Grevelingen-FM_0*_map.nc'), #zlayer
                 r'p:\1204257-dcsmzuno\2006-2012\3D-DCSM-FM\A18b_ntsu1\DFM_OUTPUT_DCSM-FM_0_5nm\DCSM-FM_0_5nm_0*_map.nc', #fullgrid
-                r'p:\dflowfm\maintenance\JIRA\05000-05999\05477\c103_ws_3d_fourier\DFM_OUTPUT_westerscheldt01_0subst\westerscheldt01_0subst_map.nc', #zsigma model without fullgrid output but with new ocean_sigma_z_coordinate variable
+                #r'p:\dflowfm\maintenance\JIRA\05000-05999\05477\c103_ws_3d_fourier\DFM_OUTPUT_westerscheldt01_0subst\westerscheldt01_0subst_map.nc', #zsigma model without fullgrid output but with new ocean_sigma_z_coordinate variable #TODO: fails since https://github.com/Deltares/xugrid/issues/68#issuecomment-1594362969 was fixed, implement/catch edge validator?
                 r'p:\archivedprojects\11206813-006-kpp2021_rmm-2d\C_Work\31_RMM_FMmodel\computations\model_setup\run_207\results\RMM_dflowfm_0*_map.nc', #2D model
                 r'p:\archivedprojects\11203379-005-mwra-updated-bem\03_model\02_final\A72_ntsu0_kzlb2\DFM_OUTPUT_MB_02\MB_02_0*_map.nc',
                 ]
@@ -167,7 +167,7 @@
     print('plot bedlevel')
     #get bedlevel and create plot with ugrid and cross section line
     fig, ax_input = plt.subplots()
-    pc = data_frommap_merged['mesh2d_flowelem_bl'].ugrid.plot(edgecolor='face',cmap='jet') #TODO: should work even better with edgecolor='none', but that results in seethrough edges anyway, report to matplotlib?
+    pc = data_frommap_merged['mesh2d_flowelem_bl'].ugrid.plot(cmap='jet') #TODO: default is edgecolor='face', should work even better with edgecolor='none', but that results in seethrough edges anyway, report to matplotlib?
     pc.set_clim(clim_bl)
     ax_input.set_aspect('equal')
     line, = ax_input.plot([], [],'o-') # empty line
@@ -195,9 +195,9 @@
         data_frommap_merged.ugrid.set_crs(crs)
         data_frommap_merged_wgs84 = data_frommap_merged.ugrid.to_crs(to_crs)
         fig, (ax1,ax2) = plt.subplots(2,1,figsize=(7,8))
-        data_frommap_merged["mesh2d_waterdepth"].isel(time=0).ugrid.plot(ax=ax1, edgecolor='face')
+        data_frommap_merged["mesh2d_waterdepth"].isel(time=0).ugrid.plot(ax=ax1)
         ctx.add_basemap(ax=ax1, source=None, crs=crs, attribution=False)
-        data_frommap_merged_wgs84["mesh2d_waterdepth"].isel(time=0).ugrid.plot(ax=ax2, edgecolor='face')
+        data_frommap_merged_wgs84["mesh2d_waterdepth"].isel(time=0).ugrid.plot(ax=ax2)
         ctx.add_basemap(ax=ax2, source=None, crs=to_crs, attribution=False)
         fig.tight_layout()
         fig.savefig(os.path.join(dir_output,f'{basename}_convertedcoords'))
@@ -206,26 +206,26 @@
     #ugrid sel via x/y
     data_frommap_merged_sel = data_frommap_merged.ugrid.sel(x=sel_slice_x,y=sel_slice_y)
     fig, ax = plt.subplots()
-    pc = data_frommap_merged_sel['mesh2d_flowelem_bl'].ugrid.plot(ax=ax, linewidth=0.5, edgecolors='face', cmap='jet')
+    pc = data_frommap_merged_sel['mesh2d_flowelem_bl'].ugrid.plot(ax=ax, linewidth=0.5, cmap='jet')
     pc.set_clim(clim_bl)
     fig.tight_layout()
     fig.savefig(os.path.join(dir_output,f'{basename}_selxyslice'))
 
 
-    if 'westerscheldt01_0subst_map' not in file_nc: #TODO: skipping for esternscheldt since contourf/contour raise "ValueError: repeats may not contain negative values." (also happens in notebook) https://github.com/Deltares/xugrid/issues/68
+    if 'westerscheldt01_0subst_map' not in file_nc: #TODO: skipping for westernscheldt since contourf/contour raise "ValueError: repeats may not contain negative values." (also happens in notebook) https://github.com/Deltares/xugrid/issues/68
         print('plot bedlevel as polycollection, contourf, contour, rasterized')
         #create fancy plots, more options at https://deltares.github.io/xugrid/examples/plotting.html
         if clim_bl is None:
             vmin = vmax = None
         else:
             vmin, vmax = clim_bl #TODO: vmin/vmax are necessary upon plot initialization (instead of pc.set_clim(clim_bl)) for proper colorbar, this is also matplotlib behaviour so not xugrid specific
         fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2,figsize=(12,7),sharex=True,sharey=True)
-        pc = data_frommap_merged['mesh2d_flowelem_bl'].ugrid.plot(ax=ax1, linewidth=0.5, edgecolors='face', cmap='jet', vmin=vmin, vmax=vmax)
+        pc = data_frommap_merged['mesh2d_flowelem_bl'].ugrid.plot(ax=ax1, linewidth=0.5, cmap='jet', vmin=vmin, vmax=vmax)
         pc = data_frommap_merged['mesh2d_flowelem_bl'].ugrid.plot.contourf(ax=ax2, levels=11, cmap='jet', vmin=vmin, vmax=vmax)
         if 'cb_3d_map' not in file_nc: #TODO: cb_3d_map fails on contour with "UserWarning: No contour levels were found within the data range." (because all bedlevels are -5m) >> colorbar gives error, is this an xugrid or matplotlib issue?
             pc = data_frommap_merged['mesh2d_flowelem_bl'].ugrid.plot.contour(ax=ax3, levels=11, cmap='jet', vmin=vmin, vmax=vmax, add_colorbar=True)
-        bl_raster = dfmt.rasterize_ugrid(data_frommap_merged[['mesh2d_flowelem_bl']],resolution=raster_res) #rasterize ugrid
-        pc = bl_raster['mesh2d_flowelem_bl'].plot(ax=ax4, cmap='jet', vmin=vmin, vmax=vmax) #plot with non-ugrid method
+        bl_raster = dfmt.rasterize_ugrid(data_frommap_merged['mesh2d_flowelem_bl'],resolution=raster_res) #rasterize ugrid uds/uda
+        pc = bl_raster.plot(ax=ax4, cmap='jet', vmin=vmin, vmax=vmax) #plot with non-ugrid method
         fig.tight_layout()
         fig.savefig(os.path.join(dir_output,f'{basename}_gridbedcontour'))
 
@@ -235,7 +235,7 @@
     data_frommap_merged['mesh2d_s1_filt'] = data_frommap_merged['mesh2d_s1'].where(~bool_drycells)
     print('plot grid and values from mapdata (waterlevel on layer, 2dim, on cell centers)')
     fig, ax = plt.subplots()
-    pc = data_frommap_merged['mesh2d_s1_filt'].isel(time=timestep).ugrid.plot(edgecolor='face',cmap='jet')
+    pc = data_frommap_merged['mesh2d_s1_filt'].isel(time=timestep).ugrid.plot(cmap='jet')
     ax.set_aspect('equal')
     fig.tight_layout()
     fig.savefig(os.path.join(dir_output,f'{basename}_mesh2d_s1_filt'))
@@ -253,7 +253,7 @@
 
     print('plot grid and values from mapdata (salinity on layer, 3dim, on cell centers), on layer')
     fig, ax = plt.subplots()
-    pc = data_frommap_merged['mesh2d_sa1'].isel(time=timestep, mesh2d_nLayers=layno, nmesh2d_layer=layno, missing_dims='ignore').ugrid.plot(edgecolor='face',cmap='jet') #missing_dims='ignore' ignores .isel() on mesh2d_nLayers/nmesh2d_layer if that dimension is not present
+    pc = data_frommap_merged['mesh2d_sa1'].isel(time=timestep, mesh2d_nLayers=layno, nmesh2d_layer=layno, missing_dims='ignore').ugrid.plot(cmap='jet') #missing_dims='ignore' ignores .isel() on mesh2d_nLayers/nmesh2d_layer if that dimension is not present
     pc.set_clim(clim_sal)
     ax.set_aspect('equal')
     fig.tight_layout()
@@ -264,7 +264,7 @@
     data_frommap_timesel = data_frommap_merged.isel(time=timestep)
     data_frommap_timesel_atdepths = dfmt.get_Dataset_atdepths(data_xr=data_frommap_timesel, depths=-4, reference='z0') #depth w.r.t. z0/waterlevel/bedlevel (also possible to provide list of floats)
     fig, ax = plt.subplots()
-    pc = data_frommap_timesel_atdepths['mesh2d_sa1'].ugrid.plot(edgecolor='face',cmap='jet') #TODO: dask\array\reductions.py:640: RuntimeWarning: All-NaN slice encountered
+    pc = data_frommap_timesel_atdepths['mesh2d_sa1'].ugrid.plot(cmap='jet') #TODO: dask\array\reductions.py:640: RuntimeWarning: All-NaN slice encountered
     pc.set_clim(clim_sal)
     ax.set_aspect('equal')
     fig.tight_layout()
@@ -297,7 +297,7 @@
         data_frommap_fou_atdepth['magn_mean'] = np.sqrt(ux_mean**2+uy_mean**2).assign_attrs(magn_mean_attrs)
 
         fig,ax = plt.subplots(figsize=(9,5))
-        pc = data_frommap_fou_atdepth['magn_mean'].ugrid.plot(edgecolor='face')
+        pc = data_frommap_fou_atdepth['magn_mean'].ugrid.plot()
         fou_raster = dfmt.rasterize_ugrid(data_frommap_fou_atdepth,resolution=raster_res) #TODO: add rasterize+quiver to notebook
         fou_raster.plot.quiver(x='mesh2d_face_x',y='mesh2d_face_y',u=fou_varname_u,v=fou_varname_v,color='w',scale=5,add_guide=False)
         pc.set_clim(0,0.10)

tests/examples_workinprogress/workinprogress_curved_quiver.py

@@ -44,7 +44,7 @@
     #plot
     fig,ax = plt.subplots()
     uds_sel_speed = np.sqrt(uds_sel.ucx*uds_sel.ucx + uds_sel.ucy*uds_sel.ucy)
-    uds_sel_speed.ugrid.plot(edgecolor='face',cmap='jet',vmax=1)
+    uds_sel_speed.ugrid.plot(cmap='jet',vmax=1)
     #ax.quiver(X, Y, U, V, color='w')
     speed = np.sqrt(U*U + V*V)
     strm = dfmt.velovect(ax, X, Y, U, V, color='w', cmap='winter', arrowstyle='fancy', linewidth=speed*2, integration_direction='forward',

tests/examples_workinprogress/workinprogress_delft3D4_netcdf.py

@@ -117,7 +117,7 @@ def preprocess_delft3d4(ds):
     ax.set_title(f't={timestep} ({timestr})')
     ax.set_aspect('equal')
     pc = ax.quiver(data_nc_XZ[::2,::2], data_nc_YZ[::2,::2], vel_x[::2,::2], vel_y[::2,::2], vel_magn[::2,::2],
-              scale=3,color='w',width=0.005, edgecolor='face', cmap='jet')
+              scale=3,color='w',width=0.005, cmap='jet')
     pc.set_clim([0,0.15])
     cbar = fig.colorbar(pc, ax=ax)
     cbar.set_label('velocity magnitude (%s)'%(ds_thiery_tsel.U1.attrs['units']))

tests/examples_workinprogress/workinprogress_exporttoshapefile.py

@@ -68,7 +68,7 @@
         newdata[varname] = data_sel_var.to_numpy() #can only have faces dimension (no time/layer)
 
         fig, ax = plt.subplots()
-        data_sel_var.ugrid.plot(cmap='viridis',edgecolor='face')
+        data_sel_var.ugrid.plot(cmap='viridis')
         fig.tight_layout()
 
     timestamp = data_map_timesel.time.dt.strftime('%Y%m%d').data

tests/examples_workinprogress/workinprogress_map_monthlymean.py

@@ -52,7 +52,7 @@
 
 print('plot grid and values from mapdata (salinity on layer, 3dim, on cell centers)')
 fig, ax = plt.subplots()
-pc = data_frommap['mesh2d_sa1'].isel(nmesh2d_layer=33,month=0,time=0,missing_dims='ignore').ugrid.plot(linewidth=0.5, cmap="jet", edgecolor='face')
+pc = data_frommap['mesh2d_sa1'].isel(nmesh2d_layer=33,month=0,time=0,missing_dims='ignore').ugrid.plot(linewidth=0.5, cmap="jet")
 pc.set_clim([28,30.2])
 ax.set_aspect('equal')
 fig.tight_layout()

tests/examples_workinprogress/workinprogress_mapfile_to_regulargrid.py

@@ -83,7 +83,7 @@
 lonvar = xr.DataArray(reg_x_vec, dims=('longitude'), attrs=lonvar_attrs)
 data_xr_out['longitude'] = lonvar
 
-latvar_attrs = {'axis':'Y',
+latvar_attrs = {'axis':'Y', #TODO: attrs are not on data_rasterized.latitude
                 'reference':'geographical coordinates, WGS84 projection',
                 'units':'degrees_north',
                 '_CoordinateAxisType':'Lat',
@@ -101,7 +101,7 @@
 # data_xr_out['depth'] = depthvar
 # data_xr_out = data_xr_out.set_coords('depth')
 
-
+data_frommap_merged = data_frommap_merged[varname_list] #speeds up rasterization significantly
 data_xr_out_temp = data_xr_out.rename({'longitude':'x','latitude':'y'}) #TODO: make rasterize_like more flexible for different xy-varnames? (look at xr.interp_like)
 data_rasterized = dfmt.rasterize_ugrid(uds=data_frommap_merged, ds_like=data_xr_out_temp) #TODO: speed up by applying uds.ugrid.sel(x/y) first
 data_rasterized = data_rasterized.rename({'x':'longitude','y':'latitude'})
@@ -115,6 +115,8 @@
     #source = ctx.providers.Esri.WorldImagery
     #ctx.add_basemap(ax, attribution=False, crs='epsg:4326', source=source)
     fig.savefig(os.path.join(dir_output, f'{model}_{varname}'))
+
+print('writing to netcdf')
 data_xr_out.to_netcdf(file_nc_out)
 
 

tests/examples_workinprogress/workinprogress_waqua_netcdf_convertedwith_getdata.py

@@ -138,7 +138,7 @@
     data_nc_VELU_thin = data_nc_VELU.loc[::thinning,::thinning]
     data_nc_VELV_thin = data_nc_VELV.loc[::thinning,::thinning]
     ax.quiver(data_nc_VELU_thin.grid_x_cen, data_nc_VELU_thin.grid_y_cen, data_nc_VELU_thin, data_nc_VELV_thin, 
-              color='w',scale=10)#,width=0.005)#, edgecolor='face', cmap='jet')
+              color='w',scale=10)#,width=0.005)#, cmap='jet')
     if RMM_name=='RMM':
         ax.set_xlim([61000, 72000])
         ax.set_ylim([438000, 446000])