Merge pull request #249 from bayonato89/mapping_GEAR1

Mapping gear1

csep/utils/readers.py

@@ -1,4 +1,4 @@
-import datetime
+import datetime, time
 import math
 import re
 import warnings
@@ -884,4 +884,227 @@ def quadtree_csv_loader(csv_fname):
  rates = rates.astype(float)
  region = QuadtreeGrid2D.from_quadkeys(quadkeys, magnitudes=mws)
 
- return rates, region, mws
+ return rates, region, mws
+
+def mapping_GEAR1(GEAR1_file, area_file, grid_file, b_value=False):
+ """ Projects and extrapolates annual estimates of M 5.95+, depth <= 70 km global seismicity, 
+ provided by the Global Earthquake Activity Rate (GEAR1) model, onto any geographical region 
+ provided its longitude and latitude coordinates and a regional mean b-value.
+ 
+ 
+ Args:
+ GEAR1_file (str): Original text file containing GEAR1 annual M 5.95+ earthquake rates.
+ 
+ area_file (str): Text file containing calculated areas (in m2) of every 0.1 x 0.1 cell 
+ on Earth. These estimates are useful to report expected number of earth-
+ quakes per unit time per unit area.
+ 
+ grid_file (str): Two-column, n-row array containing the centered longitude and latitude
+ coordinates of all the cells (with spatial resolution of 0.1 x 0.1) that 
+ make up the new testing region. E.g. 25.15, 33.25
+ 25.25, 33.25
+ 25.35, 33.25
+ Each lon lat coordinate must be a two-decimal floating number ending in 5, 
+ i.e., the cell midpoint, as seen above.
+ 
+ b_value (float): If modellers wish to extroplate GEAR1 M 5.95+ earthquake rates to a lower 
+ magnitude threshold, they must provide a generic b value of the region.
+ The default magnitude threshold for extrapolating global rates is 4.95.
+ Conversely, if modelers do not wish to extrapolate earthquake rates to 
+ lower magnitudes but only project GEAR1 over a geographic region, the 
+ b_value variable should remain False.
+
+ Returns
+ GEAR1_region.dat: Input text file containing GEAR 1 earthquake rates in cells defined within 
+ the desired geographic region. This file feeds a so-called read_GEAR1_format
+ function, which translates the format in which the GEAR1 forecasts were 
+ originally provided into a pyCSEP-friendly format.
+ 
+ area_region.dat: Input text file containing the areas of all the cells defined within the
+ geographical region. This file also feeds the read_GEAR1_format function 
+ and is aimed to express earthquake rates densities as expected number of 
+ earthquakes per year per m2.
+ 
+ 
+ """
+ print ('Reading data...')
+ bulk_dataW = np.loadtxt(GEAR1_file, skiprows=1, delimiter=',')
+ bulk_areaW = np.loadtxt(area_file, skiprows=1, delimiter=',')
+
+ # This part of the code is aimed to ensure that all lon and lat coordinates are two-digits floating 
+ # numbers. This is important, because the projection of GEAR1 onto a geographical region is basically 
+ # the intersection between two Pandas data frames.
+
+ latitudesW = []
+ longitudesW = []
+
+ for i in range(len(bulk_dataW)):
+ longitudesW.append(np.float('%.2f' % round(bulk_dataW[:,0][i],2))) 
+ latitudesW.append(np.float('%.2f' % round(bulk_dataW[:,1][i],2)))
+
+ # This is the first Pandas data frame when no extrapolations are needed:
+ if not b_value: 
+ GEAR1 = pd.DataFrame()
+ GEAR1['longitude'] = longitudesW
+ GEAR1['latitude'] = latitudesW
+ GEAR1['m595'] = bulk_dataW[:,2]
+ GEAR1['m605'] = bulk_dataW[:,3]
+ GEAR1['m615'] = bulk_dataW[:,4]
+ GEAR1['m625'] = bulk_dataW[:,5]
+ GEAR1['m635'] = bulk_dataW[:,6]
+ GEAR1['m645'] = bulk_dataW[:,7]
+ GEAR1['m655'] = bulk_dataW[:,8]
+ GEAR1['m665'] = bulk_dataW[:,9]
+ GEAR1['m675'] = bulk_dataW[:,10]
+ GEAR1['m685'] = bulk_dataW[:,11]
+ GEAR1['m695'] = bulk_dataW[:,12]
+ GEAR1['m705'] = bulk_dataW[:,13]
+ GEAR1['m715'] = bulk_dataW[:,14] 
+ GEAR1['m725'] = bulk_dataW[:,15] 
+ GEAR1['m735'] = bulk_dataW[:,16]
+ GEAR1['m745'] = bulk_dataW[:,17] 
+ GEAR1['m755'] = bulk_dataW[:,18] 
+ GEAR1['m765'] = bulk_dataW[:,19] 
+ GEAR1['m775'] = bulk_dataW[:,20] 
+ GEAR1['m785'] = bulk_dataW[:,21] 
+ GEAR1['m795'] = bulk_dataW[:,22]
+ GEAR1['m805'] = bulk_dataW[:,23]
+ GEAR1['m815'] = bulk_dataW[:,24]
+ GEAR1['m825'] = bulk_dataW[:,25]
+ GEAR1['m835'] = bulk_dataW[:,26]
+ GEAR1['m845'] = bulk_dataW[:,27] 
+ GEAR1['m855'] = bulk_dataW[:,28] 
+ GEAR1['m865'] = bulk_dataW[:,29] 
+ GEAR1['m875'] = bulk_dataW[:,30] 
+ GEAR1['m885'] = bulk_dataW[:,31] 
+ GEAR1['m895'] = bulk_dataW[:,32]
+
+ # But if extrapolations are indeed desired, this is then the first Pandas data frame:
+ else:
+ bv = b_value
+ GEAR1 = pd.DataFrame()
+ GEAR1['longitude'] = longitudesW
+ GEAR1['latitude'] = latitudesW
+ GEAR1['m495'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 4.95)))
+ GEAR1['m505'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.05)))
+ GEAR1['m515'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.15)))
+ GEAR1['m525'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.25)))
+ GEAR1['m535'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.35)))
+ GEAR1['m545'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.45)))
+ GEAR1['m555'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.55)))
+ GEAR1['m565'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.65)))
+ GEAR1['m575'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.75)))
+ GEAR1['m585'] = bulk_dataW[:,2] / (10**(-bv * (5.95 - 5.85)))
+ GEAR1['m595'] = bulk_dataW[:,2]
+ GEAR1['m605'] = bulk_dataW[:,3]
+ GEAR1['m615'] = bulk_dataW[:,4]
+ GEAR1['m625'] = bulk_dataW[:,5]
+ GEAR1['m635'] = bulk_dataW[:,6]
+ GEAR1['m645'] = bulk_dataW[:,7]
+ GEAR1['m655'] = bulk_dataW[:,8]
+ GEAR1['m665'] = bulk_dataW[:,9]
+ GEAR1['m675'] = bulk_dataW[:,10]
+ GEAR1['m685'] = bulk_dataW[:,11]
+ GEAR1['m695'] = bulk_dataW[:,12]
+ GEAR1['m705'] = bulk_dataW[:,13]
+ GEAR1['m715'] = bulk_dataW[:,14]
+ GEAR1['m725'] = bulk_dataW[:,15]
+ GEAR1['m735'] = bulk_dataW[:,16]
+ GEAR1['m745'] = bulk_dataW[:,17]
+ GEAR1['m755'] = bulk_dataW[:,18]
+ GEAR1['m765'] = bulk_dataW[:,19]
+ GEAR1['m775'] = bulk_dataW[:,20]
+ GEAR1['m785'] = bulk_dataW[:,21]
+ GEAR1['m795'] = bulk_dataW[:,22]
+ GEAR1['m805'] = bulk_dataW[:,23]
+ GEAR1['m815'] = bulk_dataW[:,24]
+ GEAR1['m825'] = bulk_dataW[:,25]
+ GEAR1['m835'] = bulk_dataW[:,26]
+ GEAR1['m845'] = bulk_dataW[:,27]
+ GEAR1['m855'] = bulk_dataW[:,28]
+ GEAR1['m865'] = bulk_dataW[:,29]
+ GEAR1['m875'] = bulk_dataW[:,30]
+ GEAR1['m885'] = bulk_dataW[:,31]
+ GEAR1['m895'] = bulk_dataW[:,32]
+
+ area = pd.DataFrame()
+ area['longitude'] = longitudesW
+ area['latitude'] = latitudesW
+ area['area'] = bulk_areaW[:,2]
+
+ # This is simply an artefact to release RAM memory for further computations:
+ bulk_dataW = []
+ latitudesW = []
+ longitudesW = []
+
+ # This is the second Pandas data frame:
+ bulk_dataR = np.loadtxt(grid_file, skiprows=0, delimiter=' ')
+
+ grid_longitudes = []
+ grid_latitudes = []
+
+ for i in range(len(bulk_dataR)):
+ grid_longitudes.append(np.float('%.2f' % round(bulk_dataR[:,0][i],2))) 
+ grid_latitudes.append(np.float('%.2f' % round(bulk_dataR[:,1][i],2))) 
+
+ polygon = pd.DataFrame()
+ polygon['longitude'] = grid_longitudes
+ polygon['latitude'] = grid_latitudes 
+
+ # And this is how the intersection between both data frames is carried out: 
+ GEAR1_region = pd.merge(polygon, GEAR1, how="inner", on=['longitude', 'latitude'])
+ area_region = pd.merge(polygon, area, how="inner", on=['longitude', 'latitude'])
+
+ GEAR1_region.to_csv('GEAR1_region.dat')
+ area_region.to_csv('areas_region.dat')
+
+ print ('The output files have been stored successfully.')
+
+
+def read_GEAR1_format(filename, area_filename, magnitudes):
+ """ Reads the format in which GEAR1 forecasts were originally created to translate them into
+ a format that is readable by pyCSEP. This function is designed to be used as a loader for 
+ the GriddedForecast.from_custom function.
+ 
+ 
+ Args:
+ filename (str): Text file containing GEAR1 seismicity rates in original format. This 
+ file is the one of the output products of the mapping_GEAR1 function.
+ 
+ area_filename (str): Text file containing calculated areas (in m2) of 0.1 lon x 0.1 lat cells.
+ This file is also an output of the mapping_GEAR1 function.
+ 
+ magnitudes (array): Array of magnitude bins in which the seismicity forecast is defined.
+ E.g., [3.95, 4.05, 4.15, ..., 8.95]
+ 
+ Returns:
+ :class:`csep.core.forecasts.GriddedForecast` 
+ """
+ t0 = time.time()
+ bulk_data = np.loadtxt(filename, skiprows=1, delimiter=',')
+
+ # Construction of the testing region:
+ lons = bulk_data[:,1]
+ lats = bulk_data[:,2]
+ coords = np.column_stack([lons, lats])
+
+ # Coordinates are given as midpoints origin should be in the 'lower left' corner:
+ r = CartesianGrid2D.from_origins(coords, magnitudes=magnitudes)
+
+ # Shape: (num_space_bins, num_mag_bins)
+ bulk_data_no_coords = bulk_data[:, 3:]
+
+ # Original GEAR1 format provides cumulative rates per meter**2
+ incremental_yrly_density = np.diff(np.fliplr(bulk_data_no_coords))
+
+ # Computing the differences, but returning array with the same shape:
+ incremental_yrly_density = np.column_stack([np.fliplr(incremental_yrly_density), bulk_data_no_coords[:,-1]])
+
+ # Read in area to denormalize back onto csep grid
+ area = np.loadtxt(area_filename, skiprows=1, delimiter=',')
+
+ # Allows to use broadcasting
+ m2_per_cell = np.reshape(area[:,-1], (len(area[:,1]), 1))
+ incremental_yrly_rates = incremental_yrly_density * m2_per_cell
+
+ return incremental_yrly_rates, r, magnitudes