From 34fa33e20bd28d75c5ca0a879b4fa9b2c788afa5 Mon Sep 17 00:00:00 2001 From: Jose Bayona Date: Thu, 1 Aug 2024 15:13:16 +0100 Subject: [PATCH] Imported a few libraries and corrected minor errors to run the mapping_GEAR1 function --- csep/utils/readers.py | 42 ++++++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/csep/utils/readers.py b/csep/utils/readers.py index 7f24f7be..0e600e57 100644 --- a/csep/utils/readers.py +++ b/csep/utils/readers.py @@ -9,6 +9,7 @@ # Third-party imports import numpy +import pandas as pd # PyCSEP imports from csep.utils.time_utils import strptime_to_utc_datetime, \ @@ -18,6 +19,7 @@ from csep.utils.iris import gcmt_search from csep.core.regions import QuadtreeGrid2D from csep.core.exceptions import CSEPIOException +from csep.core.regions import CartesianGrid2D def ndk(filename): @@ -901,11 +903,11 @@ def mapping_GEAR1(GEAR1_file, area_file, grid_file, b_value=False): 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. + make up the new testing region. E.g. 25.15 33.25 + 25.25 33.25 + 25.35 33.25 + Each lon lat coordinate, separated by single blank space, 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. @@ -915,7 +917,7 @@ def mapping_GEAR1(GEAR1_file, area_file, grid_file, b_value=False): b_value variable should remain False. Returns - GEAR1_region.dat: Input text file containing GEAR 1 earthquake rates in cells defined within + GEAR1_region.dat: Input text file containing GEAR1 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. @@ -928,8 +930,8 @@ def mapping_GEAR1(GEAR1_file, area_file, grid_file, b_value=False): """ print ('Reading data...') - bulk_dataW = np.loadtxt(GEAR1_file, skiprows=1, delimiter=',') - bulk_areaW = np.loadtxt(area_file, skiprows=1, delimiter=',') + bulk_dataW = numpy.loadtxt(GEAR1_file, skiprows=1, delimiter=',') + bulk_areaW = numpy.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 @@ -939,8 +941,8 @@ def mapping_GEAR1(GEAR1_file, area_file, grid_file, b_value=False): 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))) + longitudesW.append(float('%.2f' % round(bulk_dataW[:,0][i],2))) + latitudesW.append(float('%.2f' % round(bulk_dataW[:,1][i],2))) # This is the first Pandas data frame when no extrapolations are needed: if not b_value: @@ -1038,14 +1040,14 @@ def mapping_GEAR1(GEAR1_file, area_file, grid_file, b_value=False): longitudesW = [] # This is the second Pandas data frame: - bulk_dataR = np.loadtxt(grid_file, skiprows=0, delimiter=' ') + bulk_dataR = numpy.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))) + grid_longitudes.append(float('%.2f' % round(bulk_dataR[:,0][i],2))) + grid_latitudes.append(float('%.2f' % round(bulk_dataR[:,1][i],2))) polygon = pd.DataFrame() polygon['longitude'] = grid_longitudes @@ -1080,13 +1082,13 @@ def read_GEAR1_format(filename, area_filename, magnitudes): Returns: :class:`csep.core.forecasts.GriddedForecast` """ - t0 = time.time() - bulk_data = np.loadtxt(filename, skiprows=1, delimiter=',') + + bulk_data = numpy.loadtxt(filename, skiprows=1, delimiter=',') # Construction of the testing region: lons = bulk_data[:,1] lats = bulk_data[:,2] - coords = np.column_stack([lons, lats]) + coords = numpy.column_stack([lons, lats]) # Coordinates are given as midpoints origin should be in the 'lower left' corner: r = CartesianGrid2D.from_origins(coords, magnitudes=magnitudes) @@ -1095,16 +1097,16 @@ def read_GEAR1_format(filename, area_filename, magnitudes): 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)) + incremental_yrly_density = numpy.diff(numpy.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]]) + incremental_yrly_density = numpy.column_stack([numpy.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=',') + area = numpy.loadtxt(area_filename, skiprows=1, delimiter=',') # Allows to use broadcasting - m2_per_cell = np.reshape(area[:,-1], (len(area[:,1]), 1)) + m2_per_cell = numpy.reshape(area[:,-1], (len(area[:,1]), 1)) incremental_yrly_rates = incremental_yrly_density * m2_per_cell return incremental_yrly_rates, r, magnitudes