Skip to content

Commit

Permalink
quadtree csv reader (#186)
Browse files Browse the repository at this point in the history
* quadtree csv reader
* included depth columns in reader
* reader name change
* cartesian_grid.bounds added
* check for zero area
  • Loading branch information
khawajasim authored Feb 23, 2022
1 parent ad217ed commit d13f74f
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 6 deletions.
16 changes: 11 additions & 5 deletions csep/core/regions.py
Original file line number Diff line number Diff line change
Expand Up @@ -517,6 +517,9 @@ def __init__(self, polygons, dh, name='cartesian2d', mask=None):
# index values of polygons array into the 2d cartesian grid, based on the midpoint.
self.xs = xs
self.ys = ys
# Bounds [origin, top_right]
orgs = self.origins()
self.bounds = numpy.column_stack((orgs, orgs + dh))

def __eq__(self, other):
return self.to_dict() == other.to_dict()
Expand Down Expand Up @@ -772,13 +775,16 @@ def geographical_area_from_bounds(lon1, lat1, lon2, lat2):
Returns:
Area of cell in Km2
"""
earth_radius_km = 6371.
R2 = earth_radius_km ** 2
rad_per_deg = numpy.pi / 180.0e0
if lon1 == lon2 or lat1 == lat2:
return 0
else:
earth_radius_km = 6371.
R2 = earth_radius_km ** 2
rad_per_deg = numpy.pi / 180.0e0

strip_area_steradian = 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat1) * rad_per_deg)) \
strip_area_steradian = 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat1) * rad_per_deg)) \
- 2 * numpy.pi * (1.0e0 - numpy.cos((90.0e0 - lat2) * rad_per_deg))
area_km2 = strip_area_steradian * R2 / (360.0 / (lon2 - lon1))
area_km2 = strip_area_steradian * R2 / (360.0 / (lon2 - lon1))
return area_km2

def quadtree_grid_bounds(quadk):
Expand Down
31 changes: 30 additions & 1 deletion csep/utils/readers.py
Original file line number Diff line number Diff line change
Expand Up @@ -719,7 +719,7 @@ def _parse_datetime_to_zmap(date, time):
out['second'] = dt.second
return out

def load_quadtree_forecast(ascii_fname):
def quadtree_ascii_loader(ascii_fname):
""" Load quadtree forecasted stored as ascii text file
Note: This function is adapted form csep.forecasts.load_ascii
Expand Down Expand Up @@ -756,4 +756,33 @@ def load_quadtree_forecast(ascii_fname):
# reshape rates into correct 2d format
rates = data[:, -1].reshape(n_poly, n_mag_bins)

return rates, region, mws


def quadtree_csv_loader(csv_fname):
""" Load quadtree forecasted stored as csv file
The format expects forecast as a comma separated file, in which first column corresponds to quadtree grid cell (quadkey).
The second and thrid columns indicate depth range.
The corresponding enteries in the respective row are forecast rates corresponding to the magnitude bins.
The first line of forecast is a header, and its format is listed here:
'Quadkey', depth_min, depth_max, Mag_0, Mag_1, Mag_2, Mag_3 , ....
Quadkey is a string. Rest of the values are floats.
For the purposes of defining region objects quadkey is used.
We assume that the starting value of magnitude bins are provided in the header.
Args:
csv_fname: file name of csep forecast in csv format
Returns:
rates, region, mws (numpy.ndarray, QuadtreeRegion2D, numpy.ndarray): rates, region, and magnitude bins needed
to define QuadTree forecasts
"""

data = numpy.genfromtxt(csv_fname, dtype='str', delimiter=',')
quadkeys = data[1:, 0]
mws = data[0, 3:]
rates = data[1:, 3:]
rates = rates.astype(float)
region = QuadtreeGrid2D.from_quadkeys(quadkeys, magnitudes=mws)

return rates, region, mws

0 comments on commit d13f74f

Please sign in to comment.