gdf_to_da should return array with fill_value if shape is outside grid (

#304)

* gdf_to_da should return array with fill_value if shape is outside of grid

* cache: check length of gridintersect props

* Removed some accidentally committed verbose print statements

* Pleasing Codacy

* Please Codacy

nlmod/cache.py

@@ -438,13 +438,15 @@ def _same_function_arguments(func_args_dic, func_args_dic_cache):
             mfgrid1 = {k: v for k, v in item.mfgrid.__dict__.items() if k not in excl}
             mfgrid2 = {k: v for k, v in i2.mfgrid.__dict__.items() if k not in excl}
 
-            if not is_method_equal or mfgrid1.keys() != mfgrid2.keys():
+            is_same_length_props = all(np.all(np.size(v) == np.size(mfgrid2[k])) for k, v in mfgrid1.items())
+
+            if not is_method_equal or mfgrid1.keys() != mfgrid2.keys() or not is_same_length_props:
                 logger.info(
                     "cache was created using different gridintersect, do not use cached data"
                 )
                 return False
 
-            is_other_props_equal = all([np.all(v == mfgrid2[k]) for k, v in mfgrid1.items()])
+            is_other_props_equal = all(np.all(v == mfgrid2[k]) for k, v in mfgrid1.items())
 
             if not is_other_props_equal:
                 logger.info(

nlmod/dims/grid.py

@@ -1201,7 +1201,13 @@ def gdf_to_da(
     da : xarray DataArray
         The DataArray with the projected vector data.
     """
+    da = util.get_da_from_da_ds(ds, dims=ds.top.dims, data=fill_value)
+
     gdf_cellid = gdf_to_grid(gdf, ds, ix=ix)
+
+    if gdf_cellid.size == 0:
+        return da
+
     if min_total_overlap > 0:
         gdf_cellid["area"] = gdf_cellid.area
         area_sum = gdf_cellid[["cellid", "area"]].groupby("cellid").sum()
@@ -1232,8 +1238,6 @@ def gdf_to_da(
         else:
             gdf_agg.set_index(gdf_cellid.cellid.values, inplace=True)
 
-    da = util.get_da_from_da_ds(ds, dims=ds.top.dims, data=fill_value)
-
     if ds.gridtype == "structured":
         ixs, iys = zip(*gdf_agg.index.values)
         da.values[ixs, iys] = gdf_agg[column]
@@ -1654,7 +1658,15 @@ def gdf_to_grid(
             elif shp[geometry].geom_type in ["Polygon", "MultiPolygon"]:
                 shpn["area"] = r["areas"][i]
             shps.append(shpn)
-    gdfg = gpd.GeoDataFrame(shps, geometry=geometry, crs=gdf.crs)
+
+    if len(shps) == 0:
+        # Unable to determine the column names, because no geometries intersect with the grid
+        logger.info("No geometries intersect with the grid")
+        columns = gdf.columns.to_list() + ["area", "length", "cellid"]
+    else:
+        columns = None  # adopt from shps
+
+    gdfg = gpd.GeoDataFrame(shps, columns=columns, geometry=geometry, crs=gdf.crs)
     gdfg.index.name = gdf.index.name
     return gdfg