Filter data within a polygon

Dear all,

If I would like to plot the data within a polygon in PyGMT v0.7 under ubuntu, will be it possible?

I worked around with this code:

df = pd.read_csv('example_cat.dat', sep=',')

    polygon = '/home/asus/Documents/PyTools/Points/Sipesipe_poly_2.txt'

    # Read in polygon coordinates
    with open(polygon) as f:
        lines = f.readlines()
        coords = [[float(x) for x in line.strip().split()] for line in lines]

    # Filter data within the polygon using pandas
    mask = (df['Lon'] >= coords[0][0]) & (df['Lon'] <= coords[1][0]) & \
           (df['Lat'] >= coords[0][1]) & (df['Lat'] <= coords[2][1])
    df_within_polygon = df[mask]

I attach the polygon (Sipesipe_polyu_2.txt) and the catalog (example_cat.dat (2.9 KB) )
Sipesipe_poly_2.txt (145 Bytes)

Is there any geopandas within Pygmt to enhance the code I showed?

Thank you.
Tonino

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:

1 Like

@yvonnefroehlich, thank you for the references and code.