Repository: SCECcode/pycsep

Issue:
CartesianGrid2D not having magnitude attribute instantiation
The class representing a region do not have a explicit self.magnitude attribute. However, this attribute is called in some occasions. For example, in https://github.com/SCECcode/pycsep/blob/5f84ea97101de0439deb1e3f5c383874c7bb3801/csep/core/catalogs.py#L717

In predetermined regions, the attributed is hard-coded into the CartesianGrid2D class. 
https://github.com/SCECcode/pycsep/blob/5f84ea97101de0439deb1e3f5c383874c7bb3801/csep/core/regions.py#L58

which also should be part of the workflow of the modeler.  However, CartesianGrid2D.magnitudes should be declared at the class __init__, and have the possibility to be declared when instantiating.

Discussion:

Solution patch:
diff --git a/README.md b/README.md
index 9e086b9..929bd17 100644
--- a/README.md
+++ b/README.md
@@ -1,15 +1,17 @@
 # pyCSEP: Collaboratory for the Study of Earthquake Predictability
-![](https://i.postimg.cc/Bb60rVQP/CSEP2-Logo-CMYK.png)
-<p align=center>
-    <a target="_blank" href="https://python.org" title="Python version"><img src="https://gist.githubusercontent.com/wsavran/efce311162c32460336a4f9892218532/raw/1b9c060efd1c6e52eb53f82d4249107417d6a5ec/pycsep_python_badge.svg">
-    <a target="_blank" href="https://pypi.org/project/pycsep"><img src="https://anaconda.org/conda-forge/pycsep/badges/downloads.svg">
-    <a target="_blank" href="https://github.com/SCECcode/pycsep/actions"><img src="https://github.com/SCECcode/pycsep/actions/workflows/build-test.yml/badge.svg">
-    <a target="_blank" href="https://github.com/SCECcode/pycsep/actions"><img src="https://github.com/SCECcode/pycsep/actions/workflows/build-sphinx.yml/badge.svg">
-    <a target="_blank" href="https://codecov.io/gh/SCECcode/pycsep"><img src="https://codecov.io/gh/SCECcode/pycsep/branch/master/graph/badge.svg?token=HTMKM29MAU">
-    <a target="_blank" href="https://www.zenodo.org/badge/latestdoi/149362283"><img src="https://www.zenodo.org/badge/149362283.svg" alt="DOI"></a>
-        <a target="_blank" a style="border-width:0" href="https://doi.org/10.21105/joss.03658">
-  <img src="https://joss.theoj.org/papers/10.21105/joss.03658/status.svg" alt="DOI badge" ></a>
-</p>
+
+[![CSEP Logo](https://i.postimg.cc/Bb60rVQP/CSEP2-Logo-CMYK.png)](https://cseptesting.org)
+
+
+[![Python Version](https://gist.githubusercontent.com/wsavran/efce311162c32460336a4f9892218532/raw/1b9c060efd1c6e52eb53f82d4249107417d6a5ec/pycsep_python_badge.svg)](https://python.org)
+[![Downloads](https://anaconda.org/conda-forge/pycsep/badges/downloads.svg)](https://pypi.org/project/pycsep)
+[![Build Tests](https://github.com/SCECcode/pycsep/actions/workflows/build-test.yml/badge.svg)](https://github.com/SCECcode/pycsep/actions)
+[![Build Docs](https://github.com/SCECcode/pycsep/actions/workflows/build-sphinx.yml/badge.svg)](https://github.com/SCECcode/pycsep/actions)
+[![Coverage](https://codecov.io/gh/SCECcode/pycsep/branch/master/graph/badge.svg?token=HTMKM29MAU)](https://codecov.io/gh/SCECcode/pycsep)
+[![Zenodo](https://www.zenodo.org/badge/149362283.svg)](https://www.zenodo.org/badge/latestdoi/149362283)
+[![JOSS](https://joss.theoj.org/papers/10.21105/joss.03658/status.svg)](https://doi.org/10.21105/joss.03658)
+
+
 
 # Description:
 The pyCSEP Toolkit helps earthquake forecast model developers evaluate their forecasts with the goal of understanding
@@ -24,12 +26,13 @@ pyCSEP should:
 
 # Table of Contents:
 1. [Software Documentation](https://docs.cseptesting.org)
-2. [Installation](#installation)
-3. [Usage](#usage)
-4. [Contributing](#contributing)
-5. [Change Log](https://github.com/SCECcode/pycsep/blob/master/CHANGELOG.md)
-6. [Credits](#credits)
-7. [License](#license)
+2. [CSEP Website](https://cseptesting.org)
+3. [Installation](#installation)
+4. [Usage](#usage)
+5. [Contributing](#contributing)
+6. [Change Log](https://github.com/SCECcode/pycsep/blob/master/CHANGELOG.md)
+7. [Credits](#credits)
+8. [License](#license)
 
 # Installation:
 pyCSEP can be installed in several ways. It can be installed using conda or pip package managers or from the 
diff --git a/csep/core/regions.py b/csep/core/regions.py
index 21f7d8e..a4e7715 100644
--- a/csep/core/regions.py
+++ b/csep/core/regions.py
@@ -52,10 +52,7 @@ def california_relm_collection_region(dh_scale=1, magnitudes=None, name="relm-ca
 
     # turn points into polygons and make region object
     bboxes = compute_vertices(origins, dh)
-    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)
-
-    if magnitudes is not None:
-        relm_region.magnitudes = magnitudes
+    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, magnitudes=magnitudes, name=name)
 
     return relm_region
 
@@ -97,10 +94,8 @@ def california_relm_region(dh_scale=1, magnitudes=None, name="relm-california",
 
     # turn points into polygons and make region object
     bboxes = compute_vertices(origins, dh)
-    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)
-
-    if magnitudes is not None:
-        relm_region.magnitudes = magnitudes
+    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh,magnitudes=magnitudes,
+                                  name=name)
 
     return relm_region
 
@@ -145,10 +140,8 @@ def italy_csep_region(dh_scale=1, magnitudes=None, name="csep-italy", use_midpoi
 
     # turn points into polygons and make region object
     bboxes = compute_vertices(origins, dh)
-    italy_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)
-
-    if magnitudes is not None:
-        italy_region.magnitudes = magnitudes
+    italy_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh,
+                                   magnitudes=magnitudes, name=name)
 
     return italy_region
 
@@ -187,10 +180,8 @@ def italy_csep_collection_region(dh_scale=1, magnitudes=None, name="csep-italy-c
 
     # turn points into polygons and make region object
     bboxes = compute_vertices(origins, dh)
-    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)
-
-    if magnitudes is not None:
-        relm_region.magnitudes = magnitudes
+    relm_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, magnitudes=magnitudes,
+                                  name=name)
 
     return relm_region
 
@@ -229,10 +220,8 @@ def nz_csep_region(dh_scale=1, magnitudes=None, name="csep-nz", use_midpoint=Tru
 
     # turn points into polygons and make region object
     bboxes = compute_vertices(origins, dh)
-    nz_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)
-
-    if magnitudes is not None:
-        nz_region.magnitudes = magnitudes
+    nz_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, magnitudes=magnitudes,
+                                name=name)
 
     return nz_region
 
@@ -271,10 +260,8 @@ def nz_csep_collection_region(dh_scale=1, magnitudes=None, name="csep-nz-collect
 
     # turn points into polygons and make region object
     bboxes = compute_vertices(origins, dh)
-    nz_collection_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, name=name)
-
-    if magnitudes is not None:
-        nz_collection_region.magnitudes = magnitudes
+    nz_collection_region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh,
+                                           magnitudes=magnitudes, name=name)
 
     return nz_collection_region
 
@@ -294,9 +281,9 @@ def global_region(dh=0.1, name="global", magnitudes=None):
     lons = cleaner_range(-180.0, 179.9, dh)
     lats = cleaner_range(-90, 89.9, dh)
     coords = itertools.product(lons,lats)
-    region = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(coords, dh)], dh, name=name)
-    if magnitudes is not None:
-        region.magnitudes = magnitudes
+    region = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(coords, dh)], dh,
+                             magnitudes=magnitudes, name=name)
+
     return region
 
 def magnitude_bins(start_magnitude, end_magnitude, dmw):
@@ -585,7 +572,7 @@ class CartesianGrid2D:
     Custom regions can be easily created by using the from_polygon classmethod. This function will accept an arbitrary closed
     polygon and return a CartesianGrid class with only points inside the polygon to be valid.
     """
-    def __init__(self, polygons, dh, name='cartesian2d', mask=None):
+    def __init__(self, polygons, dh, name='cartesian2d', mask=None, magnitudes=None):
         self.polygons = polygons
         self.poly_mask = mask
         self.dh = dh
@@ -601,6 +588,7 @@ class CartesianGrid2D:
         # Bounds [origin, top_right]
         orgs = self.origins()
         self.bounds = numpy.column_stack((orgs, orgs + dh))
+        self.magnitudes = magnitudes
 
     def __eq__(self, other):
         return self.to_dict() == other.to_dict()
@@ -758,9 +746,11 @@ class CartesianGrid2D:
             dh1 = numpy.abs(lats[1]-lats[0])
             dh = numpy.max([dh1, dh2])
 
-        region = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(origins, dh)], dh, name=name)
-        if magnitudes is not None:
-            region.magnitudes = magnitudes
+        region = CartesianGrid2D(polygons=[Polygon(bbox) for bbox in compute_vertices(origins, dh)],
+                                 dh=dh,
+                                 magnitudes=magnitudes,
+                                 name=name)
+
         return region
 
     def _build_bitmask_vec(self):


Test patch:
diff --git a/tests/test_spatial.py b/tests/test_spatial.py
index 7434891..05f560a 100644
--- a/tests/test_spatial.py
+++ b/tests/test_spatial.py
@@ -57,8 +57,9 @@ class TestCartesian2D(unittest.TestCase):
         self.origins.pop(0)
         self.origins.pop(-1)
         self.num_nodes = len(self.origins)
+        self.magnitudes = numpy.array([4, 5])
         # this is kinda ugly, maybe we want to create this in a different way, class method?
-        self.cart_grid = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(self.origins, self.dh)], self.dh)
+        self.cart_grid = CartesianGrid2D([Polygon(bbox) for bbox in compute_vertices(self.origins, self.dh)], self.dh, magnitudes=self.magnitudes)
 
     def test_polygon_mask_all_masked(self):
         polygons = [Polygon(bbox) for bbox in compute_vertices(self.origins, self.dh)]