mapproject weird behaviour

Added by Romain over 3 years ago

Hi all,

I am using mapproject function to convert an ascii.xy file representing the european coastline (extracted with gshhs).
My script is pretty simple :

#0 Set params options
gmt gmtset FORMAT_FLOAT_OUT %8.8f
gmt gmtset PROJ_ELLIPSOID Clarke-1880-IGN

#1 Get the coastline points from binary gshhs file
gmt gshhg ${PATH}/gshhs_f.b -Ic -N1 > inputs/coast.dat

#2 Remove headers
sed '3d' input/coast.dat > work/coast.xy

#3 Project from WGS84 to LAMBERT93
gmt mapproject work/coast.xy -Jl3.0/46.5/44.0/49.0/1:1 -Rg -F -S -C700000/6600000 > work/coast-lamb93.xy

I have a result but mapproject seems to coarse the data as illustrated below.

Results :

head work/coast.xy

180.00000000 68.99377800
179.99677800 68.99547200
179.99516700 68.99544400
179.99322200 68.99622200
179.99161100 68.99625000
179.99013900 68.99713900
179.98663900 68.99708300
179.98519400 68.99794400

head work/coast-lamb93.xy

3423309.15840000 14816536.67840000
3423309.15840000 14816536.67840000
3423309.15840000 14816536.67840000
3423309.15840000 14816536.67840000
3423309.15840000 14816536.67840000
3423309.15840000 14816536.67840000
3423309.15840000 14815704.37120000
3423309.15840000 14815704.37120000

It seems like all decimal from the original lon/lat file – namely work/coast.xy – are not considered by the map project.

Do you guys have any guess on what's happening ?

thanks !

Romain
Ps. I tried with both 5.1.1 and 5.1.2 version from the mac bundle. Same outputs.


Replies (8)

RE: mapproject weird behaviour - Added by Paul over 3 years ago

All calculations are double precision. More decimals possible, e.g.,

echo 180.00000000 68.99377800 | gmt mapproject -Jl3.0/46.5/44.0/49.0/1:1 -Rg -F -S -C700000/6600000 --FORMAT_FLOAT_OUT=%.24g
3423309.15840000007301569    14816536.6783999986946583    

but they might just be numerical noise. Given you are giving data points 180 degrees from central meridian and you are getting "precision" back as 10^-4 meter and complaining, I guess you must have very strict requirements for this transformation!

RE: mapproject weird behaviour - Added by Romain over 3 years ago

Thanks for your reply.

With such transformation I guess I have to consider decimal up to 10^-6 or 10^-7 from the lon/lat file.

Unfortunately augmenting the decimal in FORMAT_FLOAT_OUT does not solve my issue. As you can see with attached files the option you mentioned does not fix it. My Lambert93 projected points are still coarsened than the direct output from gshhs binary file...

Is that possible to augment the precision of the calculation itself, not the representation of the output float number ?
Thanks

R.

wgs84.png - WGS84 coast projection (919 KB)

lambert93.png - Lambert93 coast projection (928 KB)

RE: mapproject weird behaviour - Added by Paul over 3 years ago

Could you post the script that makes lambert93.png? If you are plotting the same gshhg data twice (via pscoast then via gshhg) they should be the same unless your options cause some deviations. Note that the inverse -JL projection requires iterations so some precision may be lost there if that is what you are doing.

RE: mapproject weird behaviour - Added by Romain over 3 years ago

Exactly ! That's what I would like to have (pscoast equals to ones from gshhg). But I have an issue in my script for sure. Here is my script :

  1. Params
    gmt gmtset FORMAT_FLOAT_OUT %8.8f
    gmt gmtset PROJ_ELLIPSOID Clarke-1880-IGN
    gmt gmtset PROJ_LENGTH_UNIT c
  1. Projection
    gmt gshhg gshhs_f.b -Ic -N1 > coast.dat
    gmt mapproject coast.dat -Jl3.0/46.5/44.0/49.0/1:1 -Rg -F -S -C700000/6600000 > coast-lamb93.xy
  1. Plot
    gmt pscoast -Jl3.0/46.5/44.0/49.0/1:850000 -R3/5/42.25/43.60 -Df -G#d9bb7a -Cwhite -N1/0.2p,#0000FF,solid -P -K >> impact.ps
    gmt psbasemap -Jl3.0/46.5/44.0/49.0/1:850000 -R3/5/42.25/43.60 -Bf0.5a1:longitude:/f0.25a0.5:latitude:/:.title:WeSn -P -O -K >> impact.ps
    gmt psxy coast-lamb93.xy -Jx1:850000 -R700000.0000/865629.1328/6128081.8176/6279561.7280 -S+0.2c -W0.2p,black -P -O -K >> impact.ps
    gmt ps2raster -E800 -A -Tg -P impact.ps

RE: mapproject weird behaviour - Added by Paul over 3 years ago

The problem is that you are doing a 1:1 conversion on the entire planet, yielding projected distances reaching +/- 10^20 meter. There are then some cancellations that limit the precision even if double precision is used; these should not be applied when -C is set but the general code does it so that is where things go wrong.
Instead of using the global -Rg I used the more local -R and then pscoast and psxy agree perfectly. So that is your work-around.

RE: mapproject weird behaviour - Added by Paul over 3 years ago

I have created an issue #708 to address this problem.

RE: mapproject weird behaviour - Added by Paul over 3 years ago

Fixed in r14300, see discussion in #708.

RE: mapproject weird behaviour - Added by Romain over 3 years ago

Thanks ! It does work !

R.

(1-8/8)