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

Joaquim Luis <[email protected]> Fri, 1 Aug 2014 22:37:03 -0700
Newsgroups gmane.comp.gis.gmt.user
Message-ID <[email protected]>
> I can also confirm all of Kara=E2=80=99s observations. If this is a proble=
m 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=E2=80=99ll=
 have to
> try dig up an old compiled folder=E2=80=A6
>=20

This is a known issue. For those of you that need it badly, the simplest sol=
ution is to have  installed both GMT4 and GMT5. See the cookbook for how to h=
ave both running side by side.=20

Joaquim

> 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=C2=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 psx=
y
>> to make a mercator map that contains 180deg longitude and contains
>> polygons that are truncated by both the left and right map margins - I ge=
t
>> banding across the map. Depending on whether you specify your region in
>> the form 0-360 or =C2=AD180/=C2=AD180 you can get a map with polygons tha=
t 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 =C2=AD180/180
>> - 1 is digitised counterclockwise and coordinates are =C2=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 =C2=ADR region but are
>> specified either in the form =C2=AD180/180 or 0/360, and the bottom two m=
aps
>> have the same location but the =C2=ADR=C2=B9s are specified differently. I=
 have
>> annotated this.
>>=20
>> You will see that as long as the map doesn=C2=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 =C2=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=C2=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=C2=B9=
t
>> remember which version sorry =C2=AD but I can find out tonight if needed)=
 that
>> work and now they don=C2=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 t=
o
>> 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]> wrot=
e:
>>>=20
>>> Hi Kara,
>>>=20
>>>> I am using GMT 5.1.1 (r12972) on OSX Mavericks and am having a problem
>>>> with plotting polygons (I=C2=B9m using the Mercator projection). Specif=
ically
>>>> 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 =C2=B3test_0-360.pdf=C2=B2 (artefac=
t in bottom
>>>> right) with =C2=B3test_required_plotting.pdf=C2=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 th=
e
>>> 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 a=
s
>>> *.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 yo=
u
>>> search and identify a single rogue feature much much easier.
>>>=20
>>> Possible workaround: ogr2ogr from the GDAL library has a "[-clipsrc [xmi=
n
>>> 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 o=
n
>>> gmt.soest.hawaii.edu instead.
>>>=20
>>>=20
>>>>=20
>>>> To try and get around this issue I used gmtspatial and the =C2=ADC opti=
on 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=C2=B9t account for. See
>>>> "test_0-360_gmtspatial.pdf"
>>>>=20
>>>> So when gmtspatial wasn=C2=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=C2=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=C2=B9t working (a bit of trial and error with awk and bash if
>>>> statements) and then not plotting these=C5=A0 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 Sep=
t
>>>> 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 o=
n
>>> 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 ha=
ve found a bug, discuss it here first.
> To formally report bugs or request features, please register and add New I=
ssue on gmt.soest.hawaii.edu
> To unsubscribe, send the message "signoff gmt-help" to [email protected]=
ii.edu
> Note: gmt-help will become obsolete on Sept 1, 2014 - please use forum on g=
mt.soest.hawaii.edu instead.
> <new_gmt_0025.png>

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.