Plot with shapefile

Hello everyone,
I am plotting my results through following code. I want to plot it with shapefile, with high quality map. Can you please respond how I add shapefile in this code and plot high quality map.

import pygmt
import pandas as pd


earthquake_data = pd.read_csv('catalog.csv')
stationfile="stations_BB.txt"
# Set up the figure
fig = pygmt.Figure()

# Define the region of interest 
region = [89, 93, 26, 29]

# Set up the basemap with coastlines, rivers, and political boundaries
fig.basemap(region=region, projection="M6i", frame=True)
fig.coast(shorelines=True, resolution='10m', water='lightblue', land='lightgreen', borders=1)

# Plot earthquake data using circles (adjust size and color as needed)
fig.plot(
    x=earthquake_data['longitude'],
    y=earthquake_data['latitude'],
    style='c0.1c',  # Circle with a diameter of 0.1 cm
    color='red',
    label='Earthquakes'
)

# Add a legend
fig.legend(position='JMR+o0.2c', box=True)

# Add a title
fig.text(
    text="Earthquake Catalog in Bhutan Region",
    x=0.5,
    y=1.05,
    justify='CM',
    font="12p,Helvetica-Bold",
)



# Save the figure to a file or display it
fig.show()

Hi Zamir,

Welcome to the forum. Thanks for posting the PyGMT code that you are using. This example from the PyGMT Gallery can be a place to start with plotting data from shapefiles:

You need to use GeoPandas to read the shapefile format.

Hello @Zamir,

welcome to the GMT forum :slightly_smiling_face: and sorry for my delayed response!

Currently I have difficulties to understand your issue: What do you mean by “plot it with shapefile, with high quality map”? Maybe it helps if you provide the shapefile you want to use.

As already mentioned by @EJFielding you can read a shapefile by using the Python package GeoPandas, for details please see the documentation at https://geopandas.org/en/stable/docs/reference/api/geopandas.read_file.html.
The data is read into a GeoDataFrame and the next steps (regarding plotting) depend on the kind of geometry (points, lines, or polygons). There are two PyGMT gallery examples which may serve as an orientation:

@yvonnefroehlich Thanks for your response. I have plot my catalog data with simple code. I want to add the country shape file with stations to plot my catalog inside the (https://data.humdata.org/dataset/cod-ab-btn?) the shapefile. I am facing an issue how I add this shape file in the code to plot my catalog inside that.

The shapefile data available from the posted website can be plotted like this:

import pygmt
import geopandas as gpd

# Load shapefile into GeoDataFrame
# Different levels of the data are provided: adm0, adm1, adm2
# Please select for your needs
level = "adm0"
gpd_btn = gpd.read_file(
    "btn_admbnda_bnlc_20201026_SHP/btn_admbnda_" + level + "_bnlc_20201026.shp"
)

# Set up new Figure instance
fig = pygmt.Figure()

fig.coast(
    region=[89, 93, 26, 29],
    projection="M10c",
    land="lightgray",
    frame=True,
)

# Plot the polygon(s) stored in the GeoDataFrame
fig.plot(
    data=gpd_btn,
    pen="0.5p,darkgreen",
    fill="lightgreen",
)

fig.show()
# fig.savefig(fname="plot_shapefile_" + level + ".png")

Output figures:
adm0

adm1

adm2

For completness: In your posted code you set resolution="10m". To get a higher resolution you may want to use resolution="h" (high):

import pygmt

# Set up the figure
fig = pygmt.Figure()

# Set up the basemap with coastlines, rivers, and political boundaries
fig.basemap(region=[89, 93, 26, 29], projection="M6i", frame=True)
fig.coast(shorelines=True, resolution="h", water="lightblue", land="lightgreen", borders=1)

# Save the figure to a file or display it
fig.show()
# fig.savefig(fname="plot_border_higherresolution.png")

Output figure:

@yvonnefroehlich Thank you very much for your detailed response. I have download the shapefile from the link (https://data.humdata.org/dataset/cod-ab-btn?) and placed it in the same foder where I have saved the code. I have run the code but it shows the bellow error. Please how I fix this error.

self.session.start(self, **kwargs)
File “fiona/ogrext.pyx”, line 588, in fiona.ogrext.Session.start
File “fiona/ogrext.pyx”, line 143, in fiona.ogrext.gdal_open_vector
fiona.errors.DriverError: btn_admbnda_bnlc_20201026_SHP/btn_admbnda_adm0_bnlc_20201026.shp: No such file or directory

The error message says that the file cannot be found. This can have multiple reasons:

  • Did you unpack the *zip file after downloading it?
  • A typo in the file name or file path.
  • An incomplete or a wrong file path.

Yes I have extract the zip file after downloading it, I have also changed the file path. The updated code which I am using is bellow but it shows again same error.

import pygmt
import geopandas as gpd

# Load shapefile into GeoDataFrame
# Different levels of the data are provided: adm0, adm1, adm2
# Please select for your needs
level = "adm0"
shapefile_path='/home/zamir/Downloads/gpd_btn/btn_admbnda_bnlc_20201026_SHP/btn_admbnda_adm0_bnlc_20201026.shp'
gpd_btn = gpd.read_file(shapefile_path)

gpd_btn = gpd.read_file(
    f"{shapefile_path}btn_admbnda_bnlc_20201026_SHP/btn_admbnda_{level}_bnlc_20201026.shp"
)

# Set up new Figure instance
fig = pygmt.Figure()

fig.coast(
    region=[89, 93, 26, 29],
    projection="M10c",
    land="lightgray",
    frame=True,
)

# Plot the polygon(s) stored in the GeoDataFrame
fig.plot(
    data=gpd_btn,
    pen="0.5p,darkgreen",
    fill="lightgreen",
)

fig.show()
# fig.savefig(fname=f"plot_shapefile_{level}.png")

The not found error is raised because the file path is built wrongly. In your code, you repeat the part "btn_admbnda_bnlc_20201026_SHP/btn_admbnda_adm0_bnlc_20201026.shp":

level = "adm0"
shapefile_path='/home/zamir/Downloads/gpd_btn/btn_admbnda_bnlc_20201026_SHP/btn_admbnda_adm0_bnlc_20201026.shp'
gpd_btn = gpd.read_file(shapefile_path)

gpd_btn = gpd.read_file(
    f"{shapefile_path}btn_admbnda_bnlc_20201026_SHP/btn_admbnda_{level}_bnlc_20201026.shp"
)

Maybe you want:

level = "adm0"
shapefile_path='/home/zamir/Downloads/gpd_btn/'
gpd_btn = gpd.read_file(
    f"{shapefile_path}btn_admbnda_bnlc_20201026_SHP/btn_admbnda_{level}_bnlc_20201026.shp"
)

A tip for next time posting: You can format your script as code by placing three backticks in the line before and after the block with the script:

```
your script formated as code
```

@yvonnefroehlich Thank you very much for your comprehensive response and comment for posting.