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

Paul Wessel <[email protected]> Fri, 1 Aug 2014 19:57:46 -1000
Newsgroups gmane.comp.gis.gmt.user
Message-ID <[email protected]>
Kara, Sabin-

To help me debug this issue asap, please provide just ONE polygon and a =
psxy command that demonstrates the problem - that way I don't have to =
spend time trimming down the more complicated examples.  I am on a short =
leash time-wise and could need the help!

-p

On Aug 1, 2014, at 7:37 PM, Joaquim Luis <[email protected]> wrote:

>> 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
>=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 =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 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
>>>>=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.
>>>>> =
<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.

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.