I need to do some filtering of data that GMT uses for its visualizations. (I’m creating a map of political boundaries based on some geographic features.)
Is there a way to get to the data that GMT uses, but do it in numpy/pandas?
I suppose you mean the data coast uses when you speak of political boundaries. Have you looked into the gmt coast -M option?
Dump a single multisegment ASCII (or binary, see -bo ) file to standard output. No plotting occurs.
Not sure whether pygmt has the same - the manual speaks about dumping segments but I don’t see and equivalent switch (maybe some pygmt expert can comment on this). In the worst case you could mix Python and Shell or pre-process your data in Shell, then throw your Py script at it.
Cheers,
Christian
1 Like
That’s probably what I need. It’s not essential that it happen at the python level, I just need to use the data in Python.
this is BTW yet another way of displaying “nothing” for a given position: clipping out a position or a region (or clipping out everything except the position or the region) when plotting. coast is one of the modules that possesses such capability, see option -E...+c|C... and -Q, and the examples.
others are clip, mask. There is also grdclip and grdmask do different things but may also be useful for this and many other purposes. These enable modifying the grid values or creating masks based on clipping polygons or distance from clipping points.
Check e.g. examples 15, 17, also 34, Illustration Gallery — GMT 6.5.0 documentation. These are about using country polygons.
I can get this to do what I expect for political boundaries (e.g. -N1), but it doesn’t do what I’d expect for the shorelines (-W1). It only gives me some inland lakes, not the ocean/coast.
Here’s the exact command I’m using:
gmt coast -JD16.30/52.95/38/54/60 -R7/9.85/53.33/54.75 -W1 -Di -M
(Note: the dumped points are in green.)
I was trying to get a section of the coastline of Germany and Denmark. It only gives me 437 points. It doesn’t seem to matter what level I use for W.
where level is 1-4 and represent coastline, lakeshore, island-in-lake shore, and lake-in-island-in-lake shore
Curiously, if I broaden the region, I get some coastlines, but still not all the ones I would have expected…
Does the coast command have to show the entire bin that the coast points are in, in order to dump it?
If it matters, I’m trying to make a map with custom political boundaries, representing a historical time period. I need/want to merge existing physical boundaries/coastlines with other boundaries that aren’t in GMT’s lists of boundaries.
Todd,
I think you’ve misunderstood the use of the -W parameter. In most GMT plotting routines it is used to control the pen (line thickness/color/style) and doesn’t modify any data. The levels part in the -W option allows you to style rivers, lakes etc using different pen parameters when in plotting mode.
So if you modify your extraction command to something like this:
gmt coast -JD16.30/52.95/38/54/60 -R7/9.85/53.33/54.75 -Di -M -EDE,DK > test.gmt
then plotting it with
gmt plot -JD16.30/52.95/38/54/10 -R7/9.85/53.33/54.75 -W1p,blue test.gmt -png test -Bafg
gives me the following result as map:
Modifying the gmt coast -M part by using a different resolution shoreline (-D) or excluding areas based on size (-A) or also extracting rivers/lakes (-I) should get you what you want. I think you might have to check whether you need multiple calls to gmt coast when you like to extract other features (-I or -N) as I am getting error messages when combining this in one call:
coast [ERROR]: Cannot combine -E -M with -I
coast [ERROR]: Cannot combine -E -M with -N
but don’t have time to check in detail right now - getting rid of -E in this case might be the solution.
Cheers,
Christian
1 Like