pscoast -G fill missing slice in Antarctica

[btozer]

Description of the problem

Hey @PaulWessel As mentioned earlier in the week - Here is a reproduction of the fill slice missing in Antarctica.

Full script that generated the error

gmt pscoast -Rg -JG176/-41/10c -Ggrey -Wblack -Bg45 -Dl -A8000 -Xc -Yc -P

globe.ps

Full error message

A slice is not filled by -G flag fill color in Antarctica.

System information

• Operating system: MacOs Mojave 10.14.5
• Version of GMT: 5.4.5

[seisman]

Just tried the same command with gmt 6.1.0_22d70a8_2019.08.01. It's even worse.

Note there are some white lines in this plot.

[btozer]

The white lines are also in the 5.4.5 version, it's just they can't be seen at the scale of the plot.

[PaulWessel]

The missing pieces is clearly a bug, but the thin lines might be result of aliasing in gs.

[PaulWessel]

I have not figured out why it happens, but it is related to the -A8000 option. Selecting a smaller value, e.g. -A7000 works. So the workaround is to reduce the minimum area cutoff.

[PaulWessel]

Keep having to rediscover which bin is failing... For the record (in case not solving it this weekend), one bin that is failing to appear is number 641 (selected as only bin via debug -+bin option), but see -A discussion above. Here is the failing example:

gmt pscoast -Rg -JG176/-41/18c -Gred -Wblack -Bg30 -Dl -A8000 -+641 > t.ps

[PaulWessel]

Seems related to the difference between the ice line [Default] and the grounding line (-A+ag). The nodes at (-70,-80) and (-60,-80) lie between the ice and grounding lines and when we check if feature is too small (< 8000 km^2) the in-between nodes says the polygon ID is 100 (some island) and not 5 (Antarctica) and since that island is small it decides to skip this slice... Now need to figure out why the IDs pertaining to the grounding lines are used for the ice line polygon.

[PaulWessel]

We applied the fix in #1610 which seems to work fine for now. Let's see if any related issues are raised before we close this issue.

[PaulWessel]

Hi @btozer, any other issues you are finding?

[btozer]

Hey @PaulWessel, sorry I need to figure out github notifications - I didn't see you had tagged me in this message.

I had one problem with 2kml, but I think it was my mistake. It seems the -Z + options (e.g. -a for turning on/off visibility with altitude) only work for point data (I haven't tested for others) if you explicitly state a region using the -R option. If you don't bother with -R then the -Z options seem to have no effect. But I think this may be by design?

I also had a problem with sizing geovectors (-S=) using plot, but again, may be my mistake. I'll send you an e-mail with that.

[seisman]

I believe the pscoast issue was fixed in #1610. Closing.

@btozer If you still have problems with gmt2kml, please open a new issue.

[PaulWessel]

My attempt for a fix still leaves a few things to be desired. In making earth masks I found similar issues and had to relearn what happened. While it is exactly the same problem, here is perhaps a clearer description. The first plot shows a part of Antarctica at full resolution (-Df) and all features present (-A0 [Default]). Because we have two Antarctic coastlines I paint the icefront polygon in light blue and the grounding line polygon in light red. In this case everything works as expected. Notice the tiny pink islands (if the grounding line were used as coastline):

Now, if like @btozer I use -A to restrict some features, there is problems. Here I wanted to skip tiny features < 7 km^2 in area so I added -A7. The result is several tiles being left as ocean (white) in the middle of the continent:

I have framed one of the failing tiles and at each corner I have printed level [polygon ID]. The level reflects the hierarchical level of a point and goes from 0 (ocean), 1 (land), 2 (lake), 3 (island in lake), and 4 (pond in such an island). Because the NW node of the frame lies outside the grounding line but inside the iceline it is flagged with a node level of 0.5. At runtime, GMT changes those levels to either 0 or 1 depending on which Antarctica line is selected. Since the default is the icefront (-A+ai) it becomes 1 and the entire tile is on land (top figure). However, when we are skipping features we run into problems. You can see in Fig. 1 that a tiny pink island is covering the NE node (-40/-81). It has an area < 7 km^2 so it is skipped. OK, since the node level at that point used to be 1 (land) it now needs to be reset because it is now the polygon containing that island that now is containing that node. That polygon, of course, is the blue, and with us using the default coastline this is still land, so it should remain at 1 and all would be fine. That is not what happens. Here is what I (wrongly) do. For each of the four tile corners:

1. Get the ID of the polygon that contains it
2. Check if the area of the polygon (c.GSHHS_area[ID]) is < cutoff area.
3. If true then drop the node level by 1 (since the level of the polygon containing the first polygon must be one less. This change the NE node level to 0.
4. Get the ID of the parent polygon of the tiny feature. It is -1 since it is in the ocean and has no polygon around it.

Continue to determine painting color (coast) or node level (grdlandmask) and because of the wrong node levels, tiles that have no actual coastline going through them have colors determined by the node levels and that is why it ends up wrong.

Clearly, what should happen is this. For each of the four tile corners:

1. Get the ID of the polygon that contains it
2. Check if the area of the polygon (c.GSHHS_area[ID]) is < cutoff area.
3. If true then we may or may not need to make adjustments as this depends on what Antarctica coastline we have selected: a) If iceline [Default] and ID is rock-Antarctica, do nothing (it is still on land). b) If grounding line [-A+ag] subtract 1 since we now are in the ocean . In our case, this would not change the NE node level and the four corners are all 1 (land).

I have convinced myself now and will try this recipel

[PaulWessel]

Hopefully, @btozer may be able to give this a spin on some of those plots with large -A values.

[PaulWessel]

Celebrating with a cold Corona Light tonight: