Enhance mapproject to give closed domain polygons

[PaulWessel]

This PR follows up a bit on the one earlier today (#6669) by adding an optional modifier +n to -We or -Wr. This combination will build a closed polygon outline of the oblique regions (red lines). For -Wr that is the interior oblique region (left) while for -We it is the encompassing oblique region (right). I have then added a documentation figure to mapproject -W to explain the purpose of these and why they are different:

Note the solid outlines in the figures are simple geographic polygons made of just meridians and parallels so these are never created since GMT can plot these as is (e.g., plot -Ap)

Oblique polygon for an oblique projection:

gmt mapproject -R270/20/305/25+r -JOc280/25.5/22/69/5c -Wr+n > geo.txt

Polygon for a non-oblique projections encompassing polygon:

gmt mapproject -R-15/60/68/90 -JS36/90/5c -We+n > geo.txt

[joa-quim]

I'm a bit confuse still. Andreas wanted to get a -R option, which I also missed time to time for the right case, but this PR outputs a polygon only, not an -R string.

[PaulWessel]

I'm a bit confuse still. Andreas wanted to get a -R option, which I also missed time to time for the right case, but this PR outputs a polygon only, not an -R string.

If you just want a -R string you use -WE (or -WR). See examples.

[joa-quim]

OK. And the KEY? W-( is not an output key and it now can generate a polygon.

[PaulWessel]

Right, I gotta look at the KEY. Hold a sec...

[PaulWessel]

I think KEYS are fine: There is always output to stdout here >D}, and -W turns off any input as always.

[joa-quim]

Something is not write. These both come out with nothing on CLI

gmt mapproject -R-15/60/68/90 -JS36/90/5c -We+n
gmt mapproject -R270/20/305/25+r -JOc280/25.5/22/69/5c -Wr+n

and while this work on Julia the content is strange

D = mapproject(R="270/20/305/25+r", J="Oc280/25.5/22/69/5c", W="r+n10")
...
272.159205523   25.4128898844
271.426365006   23.6273106096
270.70693247    21.8222826552
270     20
GMTdataset{Float64, 2}[]

and this crashes Julia

mapproject(R="-15/60/68/90", J="S36/90/5c", W="e+n")

[PaulWessel]

Please check if on the right branch first. I cannot understand why the two CLI lines would fail since that works here. But if not on the very latest commit in that branch then I can understand it.

[joa-quim]

You really need to drop that XCode s...

[PaulWessel]

I'm running in Xcode with malloc guards on and not exceeding any limits.

[PaulWessel]

Both commands in Julia works fine here - I think you need to updated and rebuild.

[joa-quim]

I did update some 5 min ago.

[PaulWessel]

Well, all I can say is it runs on CLI and Julia here. Guard Malloc finds nothing unusual. But one of my last commits did address a memory counter issue - so hence my question.

[joa-quim]

nope, not that. git tells me that I'm on the latest change.

[PaulWessel]

OK, then I need your help for a hint for what is going on.

[joa-quim]

We have seen this before. XCode simply doesn't detect all cases like this. VS does.
It crashes in gmt_free_segment, line 8690 but this is not a big help because the true cause is on when the memory was allocated..

[joa-quim]

Interpreting this, I would say that segment->data[col] (col = 0) was not allocated with a size big enough to the data that was copied into that pointer.

[PaulWessel]

Well, I had it right the very first time, but then I decided I need to close the polygon (but was already closed) so now we got one too many, See how this goes.

[PaulWessel]

Wrong button.

[joa-quim]

Still crashes.
But if I change (mapproject line 1031)

			uint64_t dim[4] = {1, 1, 2*(Ctrl->W.nx+Ctrl->W.ny)-3, 2};

to

			uint64_t dim[4] = {1, 1, 2*(Ctrl->W.nx+Ctrl->W.ny)-2, 2};

it no longer crashes but prints the last point duplicated.

...
270.063773465   20.1663158146
270     20
270     20

[PaulWessel]

How do you get a duplicate after I remove that code?

[joa-quim]

Sorry, I had updated but maybe too quick. It not seems to work.
Let me check in Julia.

[PaulWessel]

I will make one more change just to ensure the first and last are exact. As it is the last is subject to some rounding.

[PaulWessel]

OK, I have done my last update.

[joa-quim]

So in Julia it kind of works but it's printing the result to stdout

D = mapproject(R="270/20/305/25+r", J="Oc280/25.5/22/69/5c", W="r+n10")
...
271.426365006   23.6273106096
270.70693247    21.8222826552
270     20
GMTdataset{Float64, 2}[]

The dataset, that should hold the output, is empty. The result was printed to std...

[PaulWessel]

Maybe I need to do something more here. All of mapproject writes to stdout so not sure why this is different in the wrappers. You've used mapproject for other things in Julia like this and it returnds via the memory, right?

[joa-quim]

Yes.

[PaulWessel]

How about now? I had a wrong mode passed in GMT_Write_Data and I think that is is - but I also added Ctrl->Out.file like other modules.

[joa-quim]

Yes, now it works.

[PaulWessel]

Gerat. I will go through the table produces and see if any of them do not use Out.file (without that and the parsing, any external cannot do ->OUtfile, for instance).