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.