Hello @tonino13,
probably pygmt.select (please find the documentation at pygmt.select — PyGMT) is helpful for you, whereby point 4 (“inside one of the polygons in the polygonfile”) is relevant. So far, the -F flag (please see the upstream GMT documentation at gmtselect — GMT 6.5.0 documentation) has unfortunately no alias in PyGMT. Below you find a code example based on your code and data provided above.
Code example:
import numpy as np
import pandas as pd
import pygmt
# Load data
df = pd.read_csv('example_cat.dat', sep=',')
df_resort = df[["Lon", "Lat", "id", "year", "month", "day", "hour", "minute",
"second", "Dep", "ML", "ErrOr", "ErrLat", "ErrLon", "ErrDep",
"Gap", "Np", "Ns", "Nphases"]]
# File with polygon coordinates
polygon = 'Sipesipe_poly_2.txt'
# Filter data based on polygon
# https://www.pygmt.org/dev/api/generated/pygmt.select.html
# https://docs.generic-mapping-tools.org/dev/gmtselect.html
# https://docs.generic-mapping-tools.org/dev/gmtselect.html#f
df_within_polygon = pygmt.select(
# [currently] the dataframe must contain longitude and
# latitude in the first and second column respectively
data=df_resort,
# currently there is no alias available in PyGMT for F
F=polygon,
)
# Create figure object
fig = pygmt.Figure()
# Create basic map
fig.basemap(
region=[
np.min(df['Lon'])-0.1,
np.max(df['Lon'])+0.1,
np.min(df['Lat'])-0.1,
np.max(df['Lat'])+0.1,
],
projection="M10c",
frame="afg",
)
# Plot polygon
fig.plot(
data=polygon,
pen="1p,red",
)
# Plot ALL data points
fig.plot(
x=df['Lon'],
y=df['Lat'],
style="c0.15c",
fill="lightblue",
pen="1p,black",
)
# Plot data points WITHIN polygon
fig.plot(
x=df_within_polygon['Lon'],
y=df_within_polygon['Lat'],
style="c0.15c",
fill="orange",
pen="1p,black",
)
fig.show()
# fig.savefig(fname="filter_by_polygon.png")
Output figure:
