Ensure area is computed even if no parameters are requested

mesher.py

@@ -243,6 +243,10 @@ def main():
     if use_weights and topo_weight < 0:
         raise RuntimeError("Parameter weights must equal 1")
 
+    # the empty area to trigger area computations later
+    # needs to happen here after the regularization as this doesn't have anything to regularize
+    parameter_files['area'] = {}
+
     plgs_shp = base_name + '.shp'
 
     dem = src_ds.GetRasterBand(1)
@@ -806,7 +810,7 @@ def read_config(configfile):
 
     dem_filename = X.dem_filename.strip()
     max_area = X.max_area
-    
+
     # load any user given parameter files
     parameter_files = {}
     if hasattr(X, 'parameter_files'):

pysrc/mesher_utls/MPI_do_parameterize.py

@@ -133,24 +133,6 @@ def do_parameterize(gt, is_geographic, mesh,
 
     srs_out.ImportFromProj4(srs_proj4)
 
-    # for key, data in parameter_files.items():
-    #     if data['file'] is None:
-    #         parameter_files[key]['file'] = []
-    #         for f in data['filename']:
-    #             ds = gdal.Open(f)
-    #             if ds is None:
-    #                 raise RuntimeError('Error: Unable to open raster for: %s' % key)
-    #             parameter_files[key]['file'].append(ds)
-
-    # for key, data in initial_conditions.items():
-    #     if data['file'] is None:
-    #         initial_conditions[key]['file'] = []
-    #         for f in data['filename']:
-    #             ds = gdal.Open(f)
-    #             if ds is None:
-    #                 raise RuntimeError('Error: Unable to open raster for: %s' % key)
-    #             initial_conditions[key]['file'].append(ds)
-
     params['id'] = int(elem)
 
     v0 = mesh['mesh']['elem'][elem][0]
@@ -170,6 +152,7 @@ def do_parameterize(gt, is_geographic, mesh,
     ring.AddPoint(mesh['mesh']['vertex'][v0][0],
                   mesh['mesh']['vertex'][v0][1])  # add again to complete the ring.
 
+
     # need this for the area calculation
     tpoly = ogr.Geometry(ogr.wkbPolygon)
     tpoly.AddGeometry(ring)
@@ -179,27 +162,31 @@ def do_parameterize(gt, is_geographic, mesh,
 
     mem_layer.CreateFeature(feature)
 
-    area = 0
-    # if the output is geographic, we need to project to get a reasonable area
-    if is_geographic:
-        #use an equal area mollweide projection
-        srs_moll = osr.SpatialReference()
+    # if we are computing area, we must handle this slightly differ
+    if current_parameter == 'area':
+        area = 0
+        # if the output is geographic, we need to project to get a reasonable area
+        if is_geographic:
+            #use an equal area mollweide projection
+            srs_moll = osr.SpatialReference()
 
-        srs_moll.ImportFromProj4("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ")
+            srs_moll.ImportFromProj4("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs ")
 
-        # go from what we are outputting (which is geographic) to moll to compute the area
-        transform = osr.CoordinateTransformation(srs_out, srs_moll)
-        p = tpoly.Clone()
-        p.Transform(transform)
+            # go from what we are outputting (which is geographic) to moll to compute the area
+            transform = osr.CoordinateTransformation(srs_out, srs_moll)
+            p = tpoly.Clone()
+            p.Transform(transform)
 
-        area = p.GetArea()
-    else:
-        area = tpoly.GetArea()
+            area = p.GetArea()
+        else:
+            area = tpoly.GetArea()
 
-    params['area'] = area
+        params[current_parameter] = area
+        return params  # exit early
 
     # calculate new geotransform of the feature subset
     geom = feature.geometry()
+
     src_offset = bbox_to_pixel_offsets(gt, geom.GetEnvelope(), RasterXSize, RasterYSize)
     new_gt = (
         (gt[0] + (src_offset[0] * gt[1])),
@@ -281,12 +268,13 @@ def main(pickle_file: str,
         parameter_files[key]['file'] = []
 
         # load the data
-        for f in data['filename']:
-            ds = gdal.Open(f)
-            if ds is None:
-                raise RuntimeError(f'Error: Unable to open raster for: {key}')
+        if key != 'area':
+            for f in data['filename']:
+                ds = gdal.Open(f)
+                if ds is None:
+                    raise RuntimeError(f'Error: Unable to open raster for: {key}')
 
-            parameter_files[key]['file'].append(ds)
+                parameter_files[key]['file'].append(ds)
 
         for elem in range(0, mesh['mesh']['nelem']):
             ret = do_parameterize(gt, is_geographic, mesh, parameter_files, key,

setup.py

@@ -45,7 +45,7 @@ def get_installed_gdal_version():
 USE_CONAN = str(USE_CONAN).upper() 
 
 setup(name='mesher',
-    version='2.0.2',
+    version='2.0.3',
     description='Landsurface model mesh generation',
     long_description="""
     Mesher is a novel multi-objective unstructured mesh generation software that allows mesh generation to be generated from an arbitrary number of hydrologically important features while maintaining a variable spatial resolution.