Re: Problem with plotting polygons using PSXY and also problem with GMTSPATIAL

Paul Wessel <[email protected]> Tue, 5 Aug 2014 13:54:19 -1000
Newsgroups gmane.comp.gis.gmt.user
Message-ID <[email protected]>
Hi guys-

I think I have fixed this in r13414 (GMT5) and r10251 (GMT4).  I do not =
have time to test things thoroughly today though.  Could you please give =
it a shot and if you still find clipping issue then please add a new =
issue on the tracker with a small example showing failure.

Cheers, Paul

On Aug 2, 2014, at 7:49 PM, Sabin Zahirovic =
<[email protected]> wrote:

> Hi Joaquim,
>=20
> I do have GMT4 and GMT5 running side by side. However, the bug is in =
both
> GMT4 and GMT5. The workaround is to install an older version of GMT4, =
such
> as GMT 4.5.8.=20
>=20
> Paul, thanks for offering to look at this. Your time and help is much
> appreciated!
>=20
> Cheers,
> Sabin
>=20
> --
>=20
> MR SABIN ZAHIROVIC | PhD Candidate
> School of Geosciences | Faculty of Science
>=20
>=20
>=20
>=20
> THE UNIVERSITY OF SYDNEY
> Rm 414, Madsen Building F09 | The University of Sydney | NSW | 2006
> M +61 416 775 589 =20
> E [email protected] | W http://www.earthbyte.org
> <http://www.earthbyte.org/>
>=20
> CRICOS 00026A
> This email plus any attachments to it are confidential. Any =
unauthorised
> use is strictly prohibited. If you receive this email in error, please
> delete it and any attachments.
>=20
>=20
>=20
>=20
> On 2/08/2014 3:37 pm, "Joaquim Luis" <[email protected]> wrote:
>=20
>>> I can also confirm all of Kara=92s observations. If this is a =
problem that
>>> is unlikely to be fixed in the near term, is it relatively
>>> straightforward
>>> to extract GMT 4.5.8 (or similar) from the GMT SVN? Otherwise I=92ll =
have
>>> to
>>> try dig up an old compiled folder=85
>>>=20
>>=20
>> This is a known issue. For those of you that need it badly, the =
simplest
>> solution is to have  installed both GMT4 and GMT5. See the cookbook =
for
>> how to have both running side by side.
>>=20
>> Joaquim
>>=20
>>> Thanks for the help and time in advance!
>>>=20
>>> Cheers,
>>> Sabin
>>>=20
>>> PS: I attach a sample of my own problem, the grey polygons (LIPs) =
that
>>> cross the east map boundary cause the artefact. However, as Kara
>>> demonstrated, often the west boundary is also affected.
>>>=20
>>>=20
>>> --
>>>=20
>>> MR SABIN ZAHIROVIC | PhD Candidate
>>> School of Geosciences | Faculty of Science
>>>=20
>>>=20
>>>=20
>>>=20
>>> THE UNIVERSITY OF SYDNEY
>>> Rm 414, Madsen Building F09 | The University of Sydney | NSW | 2006
>>> M +61 416 775 589
>>> E [email protected] | W http://www.earthbyte.org
>>> <http://www.earthbyte.org/>
>>>=20
>>> CRICOS 00026A
>>> This email plus any attachments to it are confidential. Any =
unauthorised
>>> use is strictly prohibited. If you receive this email in error, =
please
>>> delete it and any attachments.
>>>=20
>>>=20
>>>=20
>>>=20
>>>> On 1/08/2014 2:52 pm, "Kara Matthews" <[email protected]>
>>>> wrote:
>>>>=20
>>>> Hi Christian, and hi everyone,
>>>>=20
>>>> Thanks very much for your reply. The script I previously sent also
>>>> included my work-around in it so I could show what steps I had to =
take
>>>> to
>>>> overcome the problem. Sorry, I should have said which lines needed
>>>> commenting out. Sabin downloaded the developer version 5.1.2 =
(r13381,
>>>> 64
>>>> bit) and it didn=B9t work for him.
>>>>=20
>>>> I did some more testing and this is what I have concluded:
>>>> Using v5.1.1 (r12968, 64 bit) and 5.1.2 (r13381, 64 bit) I cannot =
use
>>>> psxy
>>>> to make a mercator map that contains 180deg longitude and contains
>>>> polygons that are truncated by both the left and right map margins =
- I
>>>> get
>>>> banding across the map. Depending on whether you specify your =
region in
>>>> the form 0-360 or =AD180/=AD180 you can get a map with polygons =
that are
>>>> truncated by one of the margins without banding, but not both.
>>>>=20
>>>> In the attached scripts I plot 8 rectangles (10 degrees wide) in =
each
>>>> map.
>>>> 4 are cut in half by the right margin and 4 are cut in half by the =
left
>>>> margin. The 4 rectangles at each margin differ as follows:
>>>> - 1 is digitised clockwise and coordinates are =AD180/180
>>>> - 1 is digitised counterclockwise and coordinates are =AD180/180
>>>> - 1 is digitised clockwise and coordinates are 0/360
>>>> - 1 is digitised counterclockwise and coordinates are 0/360
>>>>=20
>>>> I ran tests in 3 locations (three directories within
>>>> polygon_fill_problem.zip): eastern US, Pacific, Australia
>>>>=20
>>>> For each location, the top two maps have the same =ADR region but =
are
>>>> specified either in the form =AD180/180 or 0/360, and the bottom =
two maps
>>>> have the same location but the =ADR=B9s are specified differently. =
I have
>>>> annotated this.
>>>>=20
>>>> You will see that as long as the map doesn=B9t include 180 degrees
>>>> longitude, see attached gmt_poly_plotting_test_E-US.pdf (from
>>>> test1_E-US
>>>> example), then it is possible to make the map without banding, you =
just
>>>> have to choose the right =ADR format. However, if the map includes =
180
>>>> degrees longitude, see gmt_poly_plotting_test_AUS_2.pdf and
>>>> gmt_poly_plotting_test_PAC_2.pdf (from test2_PAC top two maps, and
>>>> test3_AUS bottom two maps), then it is not possible to have =
polygons
>>>> truncated on each side of the map without banding.
>>>>=20
>>>> In the attached directories contained within =
polygon_fill_problem.zip
>>>> all
>>>> the files are in place and you just have to run the scripts as they
>>>> are. I
>>>> am using 5.1.1 (r12972) 64-bit.
>>>>=20
>>>> Interestingly, I did some further testing and this problem does not
>>>> occur
>>>> using v4.5.8 (64-bit)!!
>>>>=20
>>>> I=B9m sorry to be a bother, but I had a bunch of plotting scripts =
that I
>>>> made late last year using an older developer version of GMT5 (I =
can=B9t
>>>> remember which version sorry =AD but I can find out tonight if =
needed)
>>>> that
>>>> work and now they don=B9t. Plus I know a couple of others in the =
group
>>>> having this problem.
>>>>=20
>>>>=20
>>>> Does anyone have any recommendations? To the developers, is there a
>>>> way to
>>>> change the GMT thing responsible for this back to the way it was in
>>>> v4.5.8?
>>>>=20
>>>> Thank you very very much for your help Christian and to anyone else =
who
>>>> can offer any suggestions.
>>>>=20
>>>> Cheers,
>>>> Kara
>>>>=20
>>>>=20
>>>>=20
>>>>=20
>>>>=20
>>>>=20
>>>>=20
>>>>> On 24/07/2014 7:40 am, "Christian Heine" =
<[email protected]>
>>>>> wrote:
>>>>>=20
>>>>> Hi Kara,
>>>>>=20
>>>>>> I am using GMT 5.1.1 (r12972) on OSX Mavericks and am having a
>>>>>> problem
>>>>>> with plotting polygons (I=B9m using the Mercator projection).
>>>>>> Specifically
>>>>>> I am getting banding across my map in instances when a polygon is
>>>>>> being
>>>>>> cropped by the map border. It seems that the polygon is wrapping =
the
>>>>>> wrong way around the world. Compare =B3test_0-360.pdf=B2 =
(artefact in
>>>>>> bottom
>>>>>> right) with =B3test_required_plotting.pdf=B2 (this is what I want =
the
>>>>>> map to
>>>>>> look like). A couple of my other colleagues at work are having =
the
>>>>>> same
>>>>>> problem.
>>>>>=20
>>>>> I just ran your test.sh script using Version 5.1.2_r13178 [64-bit] =
on
>>>>> OS
>>>>> 10.9.4 and get the attached result (test_0-360_chhei_crop.pdf). It =
is
>>>>> quite similar to the look of your 'desired' map. The only =
difference
>>>>> between your desired map and the one I generated using your script =
is
>>>>> the
>>>>> central part where there are white/unfilled  patches (missing
>>>>> polygons?)
>>>>> as opposed to filled space in your required map.
>>>>>=20
>>>>> Plotting COBreconstructed_${age}.00Ma.xy (2.pdf) and
>>>>> WARSreconstructed_${age}.00Ma.xy (1.pdf) independently results in =
the
>>>>> attached maps.
>>>>>=20
>>>>> For what it's worth it might be helpful to dump the reconstructed
>>>>> data as
>>>>> *.shp file and throw this into QGIS and do basic geometry =
checking: In
>>>>> QGIS there's the menu Vector -> Check Geometry Validity (btw this
>>>>> should
>>>>> also work on OGR GMT files). But I suspect that this is more =
something
>>>>> related to dateline crossings.
>>>>>=20
>>>>> Also, I ran your polygons through gmtspatial -Q+ for a handedness
>>>>> check -
>>>>> some few of your polygons are CW whereas the others are CCW =
handed.
>>>>> Not
>>>>> sure whether this causes GMT to misbehave (wouldn't assume so).
>>>>>=20
>>>>> When exporting from GPlates, I can highly recommend to use the GMT =
OGR
>>>>> format (*.gmt) as this also exports the GPlates Feature ID which =
lets
>>>>> you
>>>>> search and identify a single rogue feature much much easier.
>>>>>=20
>>>>> Possible workaround: ogr2ogr from the GDAL library has a =
"[-clipsrc
>>>>> [xmin
>>>>> ymin xmax ymax]" option which you can also use to clip the data.
>>>>>=20
>>>>> Sorry for not being able to help with specific debugging info with
>>>>> your
>>>>> setup. Maybe an update to the recent SVN versions will fix some of =
the
>>>>> problems?=20
>>>>>=20
>>>>> Greetings from NL,
>>>>> Christian
>>>>>=20
>>>>>=20
>>>>>=20
>>>>> Mailing list for GMT discussions of all kinds.  If you are not =
sure
>>>>> you
>>>>> have found a bug, discuss it here first.
>>>>> To formally report bugs or request features, please register and =
add
>>>>> New
>>>>> Issue on gmt.soest.hawaii.edu
>>>>> To unsubscribe, send the message "signoff gmt-help" to
>>>>> [email protected]
>>>>> Note: gmt-help will become obsolete on Sept 1, 2014 - please use
>>>>> forum on
>>>>> gmt.soest.hawaii.edu instead.
>>>>>=20
>>>>>=20
>>>>>>=20
>>>>>> To try and get around this issue I used gmtspatial and the =ADC =
option
>>>>>> to
>>>>>> close the polygons at the map border using my -R. This sometimes
>>>>>> works,
>>>>>> but not all the time. In some cases it seems to be inventing some
>>>>>> strange polygons that I can=B9t account for. See
>>>>>> "test_0-360_gmtspatial.pdf"
>>>>>>=20
>>>>>> So when gmtspatial wasn=B9t working I used gmtconvert to separate =
my
>>>>>> polygons out into separate files to then isolate the problem
>>>>>> polygons.
>>>>>> It seems that gmtspatial was making up some polygons with strange
>>>>>> coordinates, coordinates that weren=B9t in the original polygon =
file.
>>>>>> See
>>>>>> prob_poly.txt that coincidentally has the same longitudes as my =
map
>>>>>> region.
>>>>>>=20
>>>>>> I was able to get around my problem by working out which polygons
>>>>>> weren=B9t working (a bit of trial and error with awk and bash if
>>>>>> statements) and then not plotting these=8A but my clunky =
work-around is
>>>>>> not ideal for making lots of figures, or plotting lots of polygon
>>>>>> files.
>>>>>>=20
>>>>>> I have also attached a test script and the input files in
>>>>>> problem_with_psxy.zip. I have currently uncommented the set of =
lines
>>>>>> I
>>>>>> needed to make my ideal figure. You should be able to run this =
script
>>>>>> from the folder if needed.
>>>>>>=20
>>>>>> Thanks kindly for any help you can offer!!
>>>>>>=20
>>>>>> Kind regards,
>>>>>> Kara
>>>>>>=20
>>>>>>=20
>>>>>> ---
>>>>>>=20
>>>>>> Dr Kara J. Matthews | Postdoctoral Researcher
>>>>>> School of Geosciences | Faculty of Science
>>>>>>=20
>>>>>> THE UNIVERSITY OF SYDNEY
>>>>>> Rm 403 | Madsen Building (F09) | The University of Sydney | NSW |
>>>>>> 2006
>>>>>> E: [email protected]
>>>>>> T: +61 2 9351 3625
>>>>>> W: www.earthbyte.org
>>>>>> Mailing list for GMT discussions of all kinds. If you are not =
sure
>>>>>> you
>>>>>> have found a bug, discuss it here first. To formally report bugs =
or
>>>>>> request features, please register and add New Issue on
>>>>>> gmt.soest.hawaii.edu To unsubscribe, send the message "signoff
>>>>>> gmt-help"
>>>>>> to [email protected] Note: gmt-help will become obsolete =
on
>>>>>> Sept
>>>>>> 1, 2014 - please use forum on gmt.soest.hawaii.edu instead.
>>>>>>=20
>>>>>> =
<test_0-360.pdf><test_0-360_gmtspatial.pdf><test_required_plotting.pdf
>>>>>>> <
>>>>>> p
>>>>>> rob_poly.txt><problem_with_psxy.zip>
>>>>>=20
>>>>>=20
>>>>> Mailing list for GMT discussions of all kinds.  If you are not =
sure
>>>>> you
>>>>> have found a bug, discuss it here first.
>>>>> To formally report bugs or request features, please register and =
add
>>>>> New
>>>>> Issue on gmt.soest.hawaii.edu
>>>>> To unsubscribe, send the message "signoff gmt-help" to
>>>>> [email protected]
>>>>> Note: gmt-help will become obsolete on Sept 1, 2014 - please use
>>>>> forum on
>>>>> gmt.soest.hawaii.edu instead.
>>>>=20
>>>>=20
>>>> Mailing list for GMT discussions of all kinds.  If you are not sure =
you
>>>> have found a bug, discuss it here first.
>>>> To formally report bugs or request features, please register and =
add
>>>> New
>>>> Issue on gmt.soest.hawaii.edu
>>>> To unsubscribe, send the message "signoff gmt-help" to
>>>> [email protected]
>>>> Note: gmt-help will become obsolete on Sept 1, 2014 - please use =
forum
>>>> on
>>>> gmt.soest.hawaii.edu instead.
>>>=20
>>>=20
>>> Mailing list for GMT discussions of all kinds.  If you are not sure =
you
>>> have found a bug, discuss it here first.
>>> To formally report bugs or request features, please register and add
>>> New Issue on gmt.soest.hawaii.edu
>>> To unsubscribe, send the message "signoff gmt-help" to
>>> [email protected]
>>> Note: gmt-help will become obsolete on Sept 1, 2014 - please use =
forum
>>> on gmt.soest.hawaii.edu instead.
>>> <new_gmt_0025.png>
>>=20
>> Mailing list for GMT discussions of all kinds.  If you are not sure =
you
>> have found a bug, discuss it here first.
>> To formally report bugs or request features, please register and add =
New
>> Issue on gmt.soest.hawaii.edu
>> To unsubscribe, send the message "signoff gmt-help" to
>> [email protected]
>> Note: gmt-help will become obsolete on Sept 1, 2014 - please use =
forum on
>> gmt.soest.hawaii.edu instead.
>=20
>=20
> Mailing list for GMT discussions of all kinds.  If you are not sure =
you have found a bug, discuss it here first.
> To formally report bugs or request features, please register and add =
New Issue on gmt.soest.hawaii.edu
> To unsubscribe, send the message "signoff gmt-help" to =
[email protected]
> Note: gmt-help will become obsolete on Sept 1, 2014 - please use forum =
on gmt.soest.hawaii.edu instead.

Mailing list for GMT discussions of all kinds.  If you are not sure you have found a bug, discuss it here first.
To formally report bugs or request features, please register and add New Issue on gmt.soest.hawaii.edu
To unsubscribe, send the message "signoff gmt-help" to [email protected]
Note: gmt-help will become obsolete on Sept 1, 2014 - please use forum on gmt.soest.hawaii.edu instead.