Skip to content

Commit

Permalink
Add gallery example to showcase blockmean (#1598)
Browse files Browse the repository at this point in the history
Co-authored-by: Will Schlitzer <schlitzer90@gmail.com>
Co-authored-by: Dongdong Tian <seisman.info@gmail.com>
  • Loading branch information
3 people committed Nov 21, 2021
1 parent cda8658 commit eaa214d
Showing 1 changed file with 63 additions and 0 deletions.
63 changes: 63 additions & 0 deletions examples/gallery/histograms/blockm.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
"""
Blockmean
---------
The :meth:`pygmt.blockmean` method calculates different quantities
inside blocks/bins whose dimensions are defined via the ``spacing`` parameter.
The following examples show how to calculate the averages of the given values
inside each bin and how to report the number of points inside each bin.
"""

import pygmt

# Load sample data
data = pygmt.datasets.load_japan_quakes()
# Select only needed columns
data = data[["longitude", "latitude", "depth_km"]]

# Set the region for the plot
region = [130, 152.5, 32.5, 52.5]
# Define spacing in x and y direction (150 by 150 minute blocks)
spacing = "150m"

fig = pygmt.Figure()

# Calculate mean depth in km from all events within 150x150 minute
# bins using blockmean
df = pygmt.blockmean(data=data, region=region, spacing=spacing)
# convert to grid
grd = pygmt.xyz2grd(data=df, region=region, spacing=spacing)

fig.grdimage(
grid=grd,
region=region,
frame=["af", '+t"Mean earthquake depth inside each block"'],
cmap="batlow",
)
# plot slightly transparent landmasses on top
fig.coast(land="darkgray", transparency=40)
# plot original data points
fig.plot(
x=data.longitude, y=data.latitude, style="c0.3c", color="white", pen="1p,black"
)
fig.colorbar(frame=["x+lkm"])

fig.shift_origin(xshift="w+5c")

# Calculate number of total locations within 150x150 minute bins via
# blockmean's summary parameter
df = pygmt.blockmean(data=data, region=region, spacing=spacing, summary="n")
grd = pygmt.xyz2grd(data=df, region=region, spacing=spacing)

fig.grdimage(
grid=grd,
region=region,
frame=["af", '+t"Number of points inside each block"'],
cmap="batlow",
)
fig.coast(land="darkgray", transparency=40)
fig.plot(
x=data.longitude, y=data.latitude, style="c0.3c", color="white", pen="1p,black"
)
fig.colorbar(frame=["x+lcount"])

fig.show()

0 comments on commit eaa214d

Please sign in to comment.