Skip to content

Commit

Permalink
add box weighting
Browse files Browse the repository at this point in the history
  • Loading branch information
@xumi1993
xumi1993 committed Jul 14, 2024
1 parent e635014 commit d28deaf
Show file tree
Hide file tree
Showing 2 changed files with 170 additions and 45 deletions.
213 changes: 169 additions & 44 deletions pytomoatt/src_rec.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1030,41 +1030,49 @@ def select_by_num_rec(self, num_rec: int):
"rec_points after selection: {}".format(self._count_records())
)

def select_one_event_in_each_subgrid(self, d_deg: float, d_km: float):
"""select one event in each subgrid
def _evt_group(self, d_deg: float, d_km: float):
"""group events by grid size
:param d_deg: grid size along lat and lon in degree
:type d_deg: float
:param d_km: grid size along depth axis in km
:type d_km: float
"""

self.log.SrcReclog.info(
"src_points before selection: {}".format(self.src_points.shape[0])
)
self.log.SrcReclog.info("processing... (this may take a few minutes)")

# store index of src_points as 'src_index'
self.src_points["src_index"] = self.src_points.index

# add 'lat_group' and 'lon_group' to src_points by module d_deg
self.src_points["lat_group"] = self.src_points["evla"].apply(
lambda x: int(x / d_deg)
)
self.src_points["lon_group"] = self.src_points["evlo"].apply(
lambda x: int(x / d_deg)
)

# add 'dep_group' to src_points by module d_km
self.src_points["dep_group"] = self.src_points["evdp"].apply(
lambda x: int(x / d_km)
)

# sort src_points by 'lat_group' and 'lon_group' and 'dep_group'
self.src_points = self.src_points.sort_values(
by=["lat_group", "lon_group", "dep_group"]
)

def select_one_event_in_each_subgrid(self, d_deg: float, d_km: float):
"""select one event in each subgrid
:param d_deg: grid size along lat and lon in degree
:type d_deg: float
:param d_km: grid size along depth axis in km
:type d_km: float
"""

self.log.SrcReclog.info(
"src_points before selection: {}".format(self.src_points.shape[0])
)
# self.log.SrcReclog.info("processing... (this may take a few minutes)")

# store index of src_points as 'src_index'
# self.src_points["src_index"] = self.src_points.index

# group events by grid size
self._evt_group(d_deg, d_km)

# find all events in the same lat_group and lon_group and dep_group
# and keep only on with largest nrec
self.src_points = self.src_points.groupby(
Expand All @@ -1090,6 +1098,125 @@ def select_one_event_in_each_subgrid(self, d_deg: float, d_km: float):
# self.remove_rec_by_new_src()
self.update()

def box_weighting(self, d_deg: float, d_km: float, obj="both"):
"""Weighting sources and receivers by number in each subgrid
:param d_deg: grid size along lat and lon in degree
:type d_deg: float
:param d_km: grid size along depth axis in km, (only used when obj='src' or 'both')
:type d_km: float
:param obj: Object to be weighted, options: ``src``, ``rec`` or ``both``, defaults to ``both``
:type obj: str, optional
"""
if obj == "src":
self._box_weighting_ev(d_deg, d_km)
elif obj == "rec":
self._box_weighting_st(d_deg)
elif obj == "both":
self._box_weighting_ev(d_deg, d_km)
self._box_weighting_st(d_deg)
else:
self.log.SrcReclog.error(
"Only 'src', 'rec' or 'both' are supported for obj"
)

def _box_weighting_ev(self, d_deg: float, d_km: float):
"""Weighting sources by number of sources in each subgrid
:param d_deg: grid size along lat and lon in degree
:type d_deg: float
:param d_km: grid size along depth axis in km
:type d_km: float
"""
self.log.SrcReclog.info(
"Box weighting: d_deg={}, d_km={}".format(d_deg, d_km)
)

# group events by grid size
self._evt_group(d_deg, d_km)

# count num of sources in the same lat_group and lon_group and dep_group
self.src_points["num_sources"] = self.src_points.groupby(
["lat_group", "lon_group", "dep_group"]
)["lat_group"].transform("count")

# calculate weight for each event
self.src_points["weight"] = 1 / np.sqrt(self.src_points["num_sources"])

# drop 'lat_group' and 'lon_group' and 'dep_group'
self.src_points = self.src_points.drop(
columns=["lat_group", "lon_group", "dep_group", "num_sources"]
)

def _box_weighting_st(self, d_deg: float):
"""Weighting receivers by number of sources in each subgrid
:param d_deg: grid size along lat and lon in degree
:type d_deg: float
"""
self.log.SrcReclog.info(
"Box weighting: d_deg={}".format(d_deg)
)

# group events by grid size
self.receivers["lat_group"] = self.receivers["stla"].apply(
lambda x: int(x / d_deg)
)
self.receivers["lon_group"] = self.receivers["stlo"].apply(
lambda x: int(x / d_deg)
)

# count num of sources in the same lat_group and lon_group
self.receivers["num_receivers"] = self.receivers.groupby(
["lat_group", "lon_group"]
)["lat_group"].transform("count")

# calculate weight for each event
self.receivers["weight"] = 1 / np.sqrt(self.receivers["num_receivers"])

# assign weight to rec_points
self.rec_points["weight"] = self.rec_points.apply(
lambda x: self.receivers[
(self.receivers["staname"] == x["staname"])
]["weight"].values[0],
axis=1,
)

# assign weight to rec_points_cs
# the weight is the average of the two receivers
if not self.rec_points_cs.empty:
self.rec_points_cs["weight"] = self.rec_points_cs.apply(
lambda x: np.mean([
self.receivers[
(self.receivers["staname"] == x["staname1"])
]["weight"].values[0],
self.receivers[
(self.receivers["staname"] == x["staname2"])
]["weight"].values[0]
]),
axis=1,
)

# assign weight to rec_points_cr
# the weight is the average of the one receiver and the other source
if not self.rec_points_cr.empty:
self.rec_points_cr["weight"] = self.rec_points_cr.apply(
lambda x: np.mean([
self.receivers[
(self.receivers["staname"] == x["staname"])
]["weight"].values[0],
self.src_points[
(self.src_points["event_id"] == x["event_id2"])
]["weight"].values[0]
]),
axis=1,
)

# drop 'lat_group' and 'lon_group'
self.receivers = self.receivers.drop(
columns=["lat_group", "lon_group", "num_receivers"]
)

def count_events_per_station(self):
"""
count events per station
Expand Down Expand Up @@ -1162,21 +1289,21 @@ def _generate_cs(self, max_azi_gap, max_dist_gap):
if abs(baz_values[i] - baz_values[j]) < max_azi_gap and \
abs(dist_deg_values[i] - dist_deg_values[j]) < max_dist_gap:
data_row = {
"src_index": idx,
"rec_index1": rec_indices[i],
"staname1": stanames[i],
"stla1": stlas[i],
"stlo1": stlos[i],
"stel1": stels[i],
"rec_index2": rec_indices[j],
"staname2": stanames[j],
"stla2": stlas[j],
"stlo2": stlos[j],
"stel2": stels[j],
"phase": f"{phases[i]},cs",
"tt": tts[i] - tts[j],
"weight": (weights[i] + weights[j]) / 2,
}
"src_index": idx,
"rec_index1": rec_indices[i],
"staname1": stanames[i],
"stla1": stlas[i],
"stlo1": stlos[i],
"stel1": stels[i],
"rec_index2": rec_indices[j],
"staname2": stanames[j],
"stla2": stlas[j],
"stlo2": stlos[j],
"stel2": stels[j],
"phase": f"{phases[i]},cs",
"tt": tts[i] - tts[j],
"weight": (weights[i] + weights[j]) / 2,
}
# set src_index to index
dd_data.append(data_row)
if dd_data:
Expand Down Expand Up @@ -1291,22 +1418,20 @@ def add_noise(self, range_in_sec=0.1, mean_in_sec=0.0, shape="gaussian"):
:param shape: shape of the noise distribution probability
:type shape: str. options: gaussian or uniform
"""
if shape == "uniform":
noise = (
np.random.uniform(
low=-range_in_sec, high=range_in_sec, size=self.rec_points.shape[0]
for rec_type in [self.rec_points, self.rec_points_cs, self.rec_points_cr]:
if rec_type.empty:
continue
if rec_type in [self.rec_points_cs, self.rec_points_cr]:
range_in_sec = range_in_sec * np.sqrt(2)
if shape == "uniform":
noise = np.random.uniform(
low=mean_in_sec-range_in_sec, high=mean_in_sec+range_in_sec, size=self.rec_points.shape[0]
)
elif shape == "gaussian":
noise = np.random.normal(
loc=mean_in_sec, scale=range_in_sec, size=self.rec_points.shape[0]
)
+ mean_in_sec
)
elif shape == "gaussian":
noise = np.random.normal(
loc=mean_in_sec, scale=range_in_sec, size=self.rec_points.shape[0]
)
self.rec_points["tt"] += noise
if not self.rec_points_cs.empty:
self.rec_points_cs["tt"] += noise
if not self.rec_points_cr.empty:
self.rec_points_cr["tt"] += noise
rec_type["tt"] += noise

def rotate(self, clat:float, clon:float, angle:float):
"""Rotate sources and receivers around a center point
Expand Down
2 changes: 1 addition & 1 deletion setup.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
long_description = fh.read()


VERSION = "0.2.0"
VERSION = "0.2.1"
setup(name='PyTomoATT',
version=VERSION,
author='Mijian Xu, Masaru Nagaso',
Expand Down

0 comments on commit d28deaf

Please sign in to comment.