DOCSEARCHbyAlgolia
Submit Docs

GMT.jl Julia Wrapper for Generic Mapping Tools

repository·master·Indexed 42 days ago

https://github.com/genericmappingtools/gmt.jl

A Julia wrapper for Generic Mapping Tools (GMT) providing a monolithic interface that mirrors GMT command-line syntax and a modern, keyword-argument-based Julia interface. It enables geographic and Cartesian data manipulation, including spatial averaging with blockmedian, gridding with surface and nearneighbor, and data masking with grdmask.

Tokens
2.1K
Snippets
13
Records
13
Agent score
53%

What's inside GMT.jl

  1. How GMT.jl works: Monolithic vs. Module-based interfaces

    master

    GMT.jl provides two distinct ways to interact with the Generic Mapping Tools:

    1. Monolithic Interface (Recommended)

    This approach stays close to the GMT command-line syntax. All GMT options are passed as a single text string to the gmt() command, along with any necessary data. This is the most robust method and is updated immediately as new GMT features are released.

    2. Module-based Interface (Work in Progress)

    This modern interface allows you to call GMT modules directly as Julia functions and set options using keyword arguments. While more readable and 'Julia-native', it is currently a work in progress; some GMT supplements have not yet been ported, and certain features may not work as expected.

    Example of the modern interface: Instead of the classic syntax: gmt coast -Rd -JR12 -Bg -Dc -A10000 -Gburlywood4 -Swheat1 -P > GMT_winkel.ps

    You can use the module-based approach:

    coast(region=:global, proj=:Winkel, frame=:g, area=10000, land=:burlywood4, water=:wheat1, show=true)
    coast(region=:global, proj=:Winkel, frame=:g, area=10000, land=:burlywood4, water=:wheat1, show=true)

  3. Use a system-wide GMT installation

    master

    By default, GMT.jl uses a GMT_jll artifact to provide the GMT binary on Unix/Mac, or an MSVC binary on Windows. If you prefer to use a system-wide GMT installation, follow these steps in the Julia REPL:

    1. Set the SYSTEMWIDE_GMT environment variable to 1.
    2. Rebuild the GMT package.
    3. Restart Julia.

    Note: To make this change permanent and prevent Julia from reverting to JLL artifacts during recompiles, declare SYSTEMWIDE_GMT=1 as a permanent environment variable in your shell configuration (e.g., .bashrc).

    - ENV["SYSTEMWIDE_GMT"] = 1;
- import Pkg; Pkg.build("GMT")
- restart Julia

  4. Compare GMT classic syntax with GMT.jl modern syntax

    master

    GMT.jl allows you to replace verbose GMT command-line strings with concise Julia function calls.

    Line Plot Example: Classic: gmt psxy filename -R0/10/0/5 -JX12 -W1p -Ba -P > psfile.ps GMT.jl: plot("filename", show=true)

    Coastline Example: Classic: gmt coast -Rd -JR12 -Bg -Dc -A10000 -Gburlywood4 -Swheat1 -P > GMT_winkel.ps GMT.jl: coast(region=:global, proj=:Winkel, frame=:g, area=10000, land=:burlywood4, water=:wheat1, show=true)

    coast(region=:global, proj=:Winkel, frame=:g, area=10000, land=:burlywood4, water=:wheat1, show=true)

  5. Gridding ship track bathymetry using blockmedian and surface

    master

    This workflow demonstrates how to process sparse point data (like ship track bathymetry) into a continuous grid using a two-step process:

    1. Spatial Averaging: Use blockmedian to perform robust, median-based averaging on a specified lattice (e.g., 1 arc minute) to decimate data and reduce noise.
    2. Gridding: Use surface with a minimum curvature spline in tension algorithm to create a continuous surface from the averaged data.

    This pattern is commonly used to create global datasets from irregular sampling.

    using GMT

# 1. Read in the point data (x y z ascii file)
geo = gmt("read -Td geologists.txt");

# 2. Decimate data using median spatial averaging on a 1 arc min lattice
# -R defines the region, -I defines the grid spacing
ave = gmt("blockmedian -R158:00W/156:40W/18:00N/19:40N -I1m", geo);

# 3. Grid the data using splines in tension
# -R (defaults to input region), -I (grid spacing), -T (tension parameter)
G = gmt("surface -R -I1m -T0.2", ave);

  6. Grid satellite measurements with nearest neighbor gridding

    master

    To convert irregular satellite sensor coordinates (longitude, latitude, and value) into a regular grid, use the nearneighbor command. You must provide a region (-R), an increment (-I), and a search radius (-S).

    Example usage: gmt("nearneighbor -R<region> -I<increment> -S<radius", [lon lat sst])

    G = gmt("nearneighbor -R-12/-1/33/43 -I0.01 -S0.05", [lon lat sst]);

  7. Read data from HDF5 files using gmt("read ...")

    master

    You can use the read command with the -Tg flag to extract specific datasets from HDF5 or NetCDF files. The syntax uses a URI-like format to specify the file and the internal path to the dataset.

    Example syntax for an HDF5 file: gmt("read -Tg file_path=gd?HDF5:\"file_path\"://dataset_path")

    file = "A2016152023000.L2_LAC_SST.nc=gd?HDF5:\"A2016152023000.L2_LAC_SST.nc\"";
sst  = gmt("read -Tg " * file * "://geophysical_data/sst");
qual = gmt("read -Tg " * file * "://geophysical_data/qual_sst");
lat  = gmt("read -Tg " * file * "://navigation_data/latitude");
lon  = gmt("read -Tg " * file * "://navigation_data/longitude");

  8. Generate color tables with grd2cpt and makecpt

    master

    You can create color tables (CPTs) for plotting:

    • Use grd2cpt to create a color table from an existing grid. Use --COLOR_NAN to handle NaN values.
    • Use makecpt to create a color table based on a specific scale (e.g., -Celevation).

    Example usage: cpt_s = gmt("grd2cpt -C<color_scheme> -E -M --COLOR_NAN=255", G)

    cpt_s = gmt("grd2cpt -Cjet -E -M --COLOR_NAN=255", G);
cpt_t = gmt("makecpt -Celevation -T0/3000");

  9. Perform spatial averaging with blockmedian

    master

    The blockmedian command performs robust, median-based spatial averaging. This is useful for decimating data and handling outliers in irregular datasets.

    Common flags:

    • -R<region>: Specifies the geographic region.
    • -I<spacing>: Specifies the grid spacing (e.g., 1m for 1 arc minute).

    Input can be a GMT object containing point data.

    # Example: averaging data on a 1 arc min lattice
ave = gmt("blockmedian -R158:00W/156:40W/18:00N/19:40N -I1m", geo);

  10. Convert PostScript files to PNG or PDF using psconvert

    master

    After generating a PostScript (.ps) file using GMT commands, use psconvert to transform it into other formats like PNG or PDF.

    Common flags:

    • -Tef: Convert to Encapsulated PostScript (EPS) or similar (the example uses -Tef for high-quality output).
    • -P: Portrait orientation.
    • -A: Adjust the bounding box.
    gmt("psconvert -Tef -P -A WL_example_3.ps");

  12. Create a continuous surface with surface

    master

    The surface command grids data using a minimum curvature spline in tension algorithm. It converts discrete points into a continuous grid.

    Common flags:

    • -R<region>: Specifies the region. Using -R without arguments typically inherits the region from the input data.
    • -I<spacing>: Specifies the grid spacing.
    • -T<tension>: Sets the tension parameter for the spline.

    Returns a GMT grid object containing x, y, and z coordinates.

    # Example: gridding data using splines in tension
G = gmt("surface -R -I1m -T0.2", ave);