Creating a map with contour lines

Plotting a contour map is handled by pygmt.Figure.grdcontour.

import pygmt

# Load sample earth relief data
grid = pygmt.datasets.load_earth_relief(resolution="05m", region=[-92.5, -82.5, -3, 7])

Out:

grdcut [DEBUG]: gmt_get_filename: In: -92.5/-82.5/-3/7 Out: -92.5/-82.5/-3/7
grdcut [DEBUG]: Look for file -92.5/-82.5/-3/7 in /vercel/.gmt
grdcut [DEBUG]: Look for file -92.5/-82.5/-3/7 in /vercel/.gmt/cache
grdcut [DEBUG]: Look for file -92.5/-82.5/-3/7 in /vercel/.gmt/server
grdcut [DEBUG]: Got regular w/e/s/n for region (-92.5/-82.5/-3/7)
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: Replace file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 with /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [INFORMATION]: Processing input grid
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: Object ID 0 : Registered Grid File /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 as an Input resource with geometry Surface [n_objects = 1]
grdcut [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdcut [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 1
grdcut [DEBUG]: Object ID 1 : Registered Grid Memory Copy 56149b234eb0 as an Output resource with geometry Surface [n_objects = 2]
grdcut [DEBUG]: Successfully created a new Grid container
grdcut [DEBUG]: VirtualFile name created: @GMTAPI@-S-O-G-G-G-Y-000001
grdcut [DEBUG]: GMT now running in modern mode [Session ID = 4724]
grdcut [DEBUG]: Revised options: @earth_relief_05m_p/ -R-92.5/-82.5/-3/7 -I05m -rp -G@GMTAPI@-S-O-G-G-G-Y-000001 -fg -Co+n
grdcut (gmtlib_free_tmp_arrays): tried to free unallocated memory
grdblend [DEBUG]: History: Process -R-92.5/-82.5/-3/7
grdblend [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
grdblend [DEBUG]: gmt_get_filename: In: -92.5/-82.5/-3/7 Out: -92.5/-82.5/-3/7
grdblend [DEBUG]: Look for file -92.5/-82.5/-3/7 in /vercel/.gmt
grdblend [DEBUG]: Look for file -92.5/-82.5/-3/7 in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file -92.5/-82.5/-3/7 in /vercel/.gmt/server
grdblend [DEBUG]: Got regular w/e/s/n for region (-92.5/-82.5/-3/7)
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Replace file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 with /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Map distance calculation will be using great circle approximation with authalic auxiliary latitudes and authalic (R_2) radius = 6371007.1809 m, in meter.
grdblend [DEBUG]: gmtapi_init_import: Passed family = Data Table and geometry = Text
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Object ID 2 : Registered Data Table File /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 as an Input resource with geometry Text [n_objects = 3]
grdblend (gmtapi_init_import): tried to free unallocated memory
grdblend [DEBUG]: gmtapi_init_import: Added 1 new sources
grdblend [DEBUG]: GMT_Init_IO: Returned first Input object ID = 2
grdblend [DEBUG]: GMT_Begin_IO: Mode value 1 not considered (ignored)
grdblend [DEBUG]: GMT_Begin_IO: Initialize record-by-record access for Input
grdblend [DEBUG]: gmtapi_next_io_source: Selected object 2
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000 Out: /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Given full path to file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: Found readable file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [INFORMATION]: Reading Data Table from file /vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdblend [DEBUG]: GMT_Begin_IO: Input resource access is now enabled [record-by-record]
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [DEBUG]: Chosen boundary condition for all edges: geographic
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [DEBUG]: Object ID 3 : Registered Grid Memory Reference 56149add4e90 as an Input resource with geometry Surface [n_objects = 4]
grdblend [DEBUG]: Successfully created a new Grid container
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Remote file @S90W180.earth_relief_05m_p.nc exists locally as /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Replace file S90W180.earth_relief_05m_p.nc with path /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Object ID 4 : Registered Grid File S90W180.earth_relief_05m_p.nc as an Input resource with geometry Surface [n_objects = 5]
grdblend [DEBUG]: gmtapi_import_grid: Passed ID = 4 and mode = 33
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc Out: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Call gmtgrdio_doctor_geo_increments on a geographic grid
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [DEBUG]: Grid @S90W180.earth_relief_05m_p.nc: out: -1050/1109/1163/-996 in: -1051/1110/1164/-997 skip: 0 offset: -996
grdblend [INFORMATION]: Blend file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc in -180/0/-90/90 with normal weight 1 [-996-1163]
grdblend (gmtlib_close_grd): tried to free unallocated memory
grdblend [DEBUG]: GMT_Destroy_Data: freed memory for a Grid for object 4
grdblend [DEBUG]: gmtlib_unregister_io: Unregistering object no 4 [n_objects = 4]
grdblend [DEBUG]: gmtlib_unregister_io: Object no 4 has non-NULL resource pointer
grdblend [DEBUG]: GMT_End_IO: Input resource access is now disabled
grdblend [INFORMATION]: Only 1 grid found; no blending will take place
grdblend [INFORMATION]: Processing input grids
grdblend [DEBUG]: Successfully added data array to previously registered Grid container
grdblend [DEBUG]: Remote file @S90W180.earth_relief_05m_p.nc exists locally as /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Replace file S90W180.earth_relief_05m_p.nc with path /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: gmt_get_filename: In: S90W180.earth_relief_05m_p.nc Out: S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Object ID 5 : Registered Grid File S90W180.earth_relief_05m_p.nc as an Input resource with geometry Surface [n_objects = 5]
grdblend [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdblend [DEBUG]: gmtapi_import_grid: Passed ID = 5 and mode = 33
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/cache
grdblend [DEBUG]: Look for file S90W180.earth_relief_05m_p.nc in /vercel/.gmt/server
grdblend [DEBUG]: Found readable file /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: gmt_get_filename: In: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc Out: /vercel/.gmt/server/earth/earth_relief/earth_relief_05m_p/S90W180.earth_relief_05m_p.nc
grdblend [DEBUG]: Call gmtgrdio_doctor_geo_increments on a geographic grid
grdblend [DEBUG]: Geographic input grid, longitudes span less than 360
grdblend [DEBUG]: GMT_End_IO: Input resource access is now disabled
grdblend [INFORMATION]: Processed row       0 of 120
grdblend [INFORMATION]: Processed row      10 of 120
grdblend [INFORMATION]: Processed row      20 of 120
grdblend [INFORMATION]: Processed row      30 of 120
grdblend [INFORMATION]: Processed row      40 of 120
grdblend [INFORMATION]: Processed row      50 of 120
grdblend [INFORMATION]: Processed row      60 of 120
grdblend [INFORMATION]: Processed row      70 of 120
grdblend [INFORMATION]: Processed row      80 of 120
grdblend [INFORMATION]: Processed row      90 of 120
grdblend [INFORMATION]: Processed row     100 of 120
grdblend [INFORMATION]: Processed row     110 of 120
grdblend [INFORMATION]: Processed row     120
grdblend [DEBUG]: GMT_Write_Data: Writing Grid to memory object 1 from object 3 which transfers ownership
grdblend [DEBUG]: gmtapi_begin_io: Output resource access is now enabled [container]
grdblend [DEBUG]: gmtapi_export_data: Messenger dummy output container for object 1 [item 1] freed and set resource=data=NULL
grdblend [DEBUG]: gmtapi_export_grid: Passed ID = 1 and mode = 0
grdblend [INFORMATION]: Referencing grid data to GMT_GRID memory location
grdblend [DEBUG]: Chosen boundary condition for all edges: geographic
grdblend [INFORMATION]: Set boundary condition for all edges: natural
grdblend [INFORMATION]: Set boundary condition for left   edge: natural
grdblend [INFORMATION]: Set boundary condition for right  edge: natural
grdblend [INFORMATION]: Set boundary condition for bottom edge: natural
grdblend [INFORMATION]: Set boundary condition for top    edge: natural
grdblend [DEBUG]: GMT_End_IO: Output resource access is now disabled
grdblend (gmtlib_close_grd): tried to free unallocated memory
grdblend [DEBUG]: GMT_Destroy_Data: freed memory for a Grid for object 5
grdblend [DEBUG]: gmtlib_unregister_io: Unregistering object no 5 [n_objects = 4]
grdblend [DEBUG]: gmtlib_unregister_io: Object no 5 has non-NULL resource pointer
grdblend [INFORMATION]: Blended grid size of @GMTAPI@-S-O-G-G-G-Y-000001 is 120 x 120
grdblend [INFORMATION]: All nodes assigned values
grdblend [DEBUG]: gmtlib_unregister_io: Unregistering object no 2 [n_objects = 3]
grdblend [DEBUG]: gmtlib_unregister_io: Unregistering object no 3 [n_objects = 2]
grdblend (gmtlib_free_tmp_arrays): tried to free unallocated memory
grdcut [INFORMATION]: Set boundary condition for all edges: natural
grdcut [INFORMATION]: Set boundary condition for left   edge: natural
grdcut [INFORMATION]: Set boundary condition for right  edge: natural
grdcut [INFORMATION]: Set boundary condition for bottom edge: natural
grdcut [INFORMATION]: Set boundary condition for top    edge: natural
grdcut [DEBUG]: GMT_End_IO: Input resource access is now disabled
grdcut [DEBUG]: gmtapi_begin_io: Input resource access is now enabled [container]
grdcut [DEBUG]: gmtapi_import_grid: Passed ID = 0 and mode = 2
grdcut [DEBUG]: GMT_End_IO: Input resource access is now disabled
grdcut [INFORMATION]: File spec:        W E S N dx dy n_columns n_rows:
grdcut [INFORMATION]: Old:
grdcut [INFORMATION]:   -92.5   -82.5   -3      7       0.0833333333333 0.0833333333333 120     120
grdcut [INFORMATION]: New:
grdcut [INFORMATION]:   -92.5   -82.5   -3      7       0.0833333333333 0.0833333333333 120     120
grdcut [DEBUG]: Object ID 6 : Registered Grid File /tmp/pygmt-oal3rlpl.nc as an Output resource with geometry Surface [n_objects = 3]
grdcut [DEBUG]: gmtapi_begin_io: Output resource access is now enabled [container]
grdcut [DEBUG]: gmtapi_export_grid: Passed ID = 6 and mode = 0
grdcut [INFORMATION]: Writing grid to file /tmp/pygmt-oal3rlpl.nc
grdcut [DEBUG]: region: 0 0, grid: -92.5 -82.5
grdcut [DEBUG]: -> region: -92.5 -82.5, grid: -92.5 -82.5
grdcut [DEBUG]: row: 0 119, col: 0 119
grdcut [DEBUG]: netCDF Library version: 4
grdcut [DEBUG]: packed z-range: [-4490.5,1140]
grdcut [DEBUG]: GMT_End_IO: Output resource access is now disabled
grdcut [DEBUG]: gmtlib_garbage_collection: Destroying object: C=0 A=0 ID=0 W=Input F=Grid M=File S=Used P=56149add4e90 N=/vercel/.gmt/sessions/gmt_session.4724/=tiled_84_GX.000000
grdcut [DEBUG]: GMTAPI_Garbage_Collection freed 1 memory objects
grdcut [DEBUG]: gmtlib_unregister_io: Unregistering object no 0 [n_objects = 2]
grdcut [DEBUG]: gmtlib_unregister_io: Unregistering object no 1 [n_objects = 1]
grdcut [DEBUG]: gmtlib_unregister_io: Unregistering object no 6 [n_objects = 0]
grdcut (gmtlib_free_tmp_arrays): tried to free unallocated memory
@earth_relief_05m -G/tmp/pygmt-oal3rlpl.nc -R-92.5/-82.5/-3/7 -Vd /tmp/pygmt-oal3rlpl.nc 60352

Create contour plot

The pygmt.Figure.grdcontour method takes the grid input. It plots annotated contour lines, which are thicker and have the elevation/depth written on them, and unannotated contour lines. In the example below, the default contour line intervals are 500 meters, with an annotated contour line every 1000 meters. By default, it plots the map with the equidistant cylindrical projection and with no frame.

fig = pygmt.Figure()
fig.grdcontour(grid=grid)
fig.show()

Out:

<IPython.core.display.Image object>

Contour line settings

Use the annotation and interval arguments to adjust contour line intervals. In the example below, there are contour intervals every 250 meters and annotated contour lines every 1,000 meters.

fig = pygmt.Figure()
fig.grdcontour(
    annotation=1000,
    interval=250,
    grid=grid,
)
fig.show()

Out:

<IPython.core.display.Image object>

Contour limits

The limit argument sets the minimum and maximum values for the contour lines. The argument takes the low and high values, and is either a list (as below) or a string limit="-4000/-2000".

fig = pygmt.Figure()
fig.grdcontour(
    annotation=1000,
    interval=250,
    grid=grid,
    limit=[-4000, -2000],
)
fig.show()

Out:

<IPython.core.display.Image object>

Map settings

The pygmt.Figure.grdcontour method accepts additional arguments, including setting the projection and frame.

fig = pygmt.Figure()
fig.grdcontour(
    annotation=1000,
    interval=250,
    grid=grid,
    limit=[-4000, -2000],
    projection="M10c",
    frame=True,
)
fig.show()

Out:

<IPython.core.display.Image object>

Adding a colormap

The pygmt.Figure.grdimage method can be used to add a colormap to the contour map. It must be called prior to pygmt.Figure.grdcontour to keep the contour lines visible on the final map. If the projection argument is specified in the pygmt.Figure.grdimage method, it does not need to be repeated in the pygmt.Figure.grdcontour method.

fig = pygmt.Figure()
fig.grdimage(
    grid=grid,
    cmap="haxby",
    projection="M10c",
    frame=True,
)
fig.grdcontour(
    annotation=1000,
    interval=250,
    grid=grid,
    limit=[-4000, -2000],
)
fig.show()

Out:

<IPython.core.display.Image object>