Using Julia to plot a GeoJSON coordinate matrix (GMT.jl vs gmt command)

Lounge
1

Hi,

I’ve been exploring GMT.jl again lately (thanks to Joaquim for this amazing package and his responses to my questions on github) and I thought I would put this here for reference, since I didn’t find so many examples out there. In the following Julia code I will show two methods for plotting a memory-resident data matrix (in this case just coordinates) on a global basemap. I will use GeoJSON formatted tectonic plate boundary data saved in the file plate_boundaries.json (in the same directory as the Julia script).

The first method uses Julia’s Cmd syntax to run the external gmt binary. This is unsurprisingly much slower than using GMT.jl (second method). The trick with GMT.jl is to use the mat2ds function to convert the matrix to a GMTdataset which is an allowed object type for the GMT.jl functions. It was much simpler than I had thought – I was overcomplicating things a lot until I found this function. Note that there is an even simpler method on the latest GMT.jl master using gmtread in this case, see comment below. Here goes:

using BenchmarkTools  # https://github.com/JuliaCI/BenchmarkTools.jl
using DelimitedFiles

using GMT
using GeoJSON  # https://github.com/JuliaGeo/GeoJSON.jl


function read()
    # Although GMT.jl can "read" GeoJSON via `gmtread("plate_boundaries.json", data = true)`,
    # the JSON string is not parsed into native data types.
    # I didn't really want to write a GeoJSON parser, and luckily it's already available!
    boundaries = open("plate_boundaries.json") do file
        geo2dict(GeoJSON.read(file))
    end
    return boundaries
end


function gmtcall()
    run(`gmt begin plate_boundaries`)
    run(`gmt psbasemap -Rg -JR14c -Baf -BWSen`)

    boundaries = read()

    for data in boundaries["features"]
        open(`gmt plot -Wthinner,black`, write = true) do io
            writedlm(io, reduce(hcat, data["geometry"]["coordinates"]) |> transpose |> collect)
        end
    end

    run(`gmt end`)
end


function gmtjl()
    boundaries = read()

    basemap(region = :g, projection = :Winkel)
    for data in boundaries["features"]
        plot!(
            mat2ds(reduce(hcat, data["geometry"]["coordinates"]) |> transpose |> collect),
            pen = (:thinner, :black),
        )
    end
end


function main()
    # Benchmarking in global scope is discouraged, so we do it in here.
    @btime gmtcall()  # ~ 9.6s and 21 MiB
    @btime gmtjl()  # ~ 0.4s and 20 MiB
end

main()

If you can see any obvious flaws with this code, please don’t hesitate to point them out! I’m quite new to both Juila and GMT…

Well, that’s all for now. Thanks to everyone who is working on GMT software :slight_smile:

2

@adigitoleo Thanks for those examples on how to manipulate data read by the GeoJSON package. I might need them one day … but not today :slight_smile:

One comment about your first example. If running a plain hard-core GMT syntax is wished then better use the monolithic syntax. I.e.

basemap("-Rg -JR14c -Baf -BWSen")

instead of

run(`gmt psbasemap -Rg -JR14c -Baf -BWSen`)

but we don’t need nothing of this. Like I mentioned in the issue, there must be something strange in your GMT.jl installation (too old?) because .json files are detected automatically and if not the keyword to read them is ogr=true and not data=true.

I found one plate boundaries file in json here and it works beautifully

D = gmtread("PB2002_boundaries.json");
imshow(D, proj=:guess, lc=:blue, coast=true)

GMTjl_tmp
GMTjl_tmp1817×1039 181 KB

1 Like
3

Aha, that’s even simpler! I’ll investigate my gmt more (it’s version 6.1.1).

The run() stuff is quite silly, agreed, but I was just satisfying my own curiosity to be honest.

I guess the example could still be useful in case the data is coming from Julia memory in the first place.

4

Yes, that’s what in mind with my first sentence.