GenericMappingTools / GMT.jl

Generic Mapping Tools Library Wrapper for Julia
Other
193 stars 28 forks source link

grdimage incorrectly handling longitude coordinates #1493

Closed melodyjulia closed 1 month ago

melodyjulia commented 1 month ago

I tried the following code snippet with the latest GMT.jl version 1.16.0, it shows wrong latitude lines:

lon = 0:1:359
lat = -89.5:1:89.5
var = zeros(360,180)
g = mat2grid(Array(var'), x = lon, y = lat)

c = makecpt(color=:rainbow, range=(-1,1,0.2))

grdimage(g, proj=:Winkel, colorbar=false, coast=true, frame=(annot="60d", grid="60d", title="Winkel"),
         color=c, par=(MAP_FRAME_TYPE=:plain,FONT=4,FONT_TITLE=8))

showfig()

If the longitude coordinates are changed to the following:

lon = 0.5:1:359.5

the latitude lines look reasonable.

image
joa-quim commented 1 month ago

If the longitude coordinates are changed to the following:

lon = 0.5:1:359.5

Why you say that?

lon = 0:1:359; lat = -89.5:1:89.5; var = zeros(Float32,180,360);

G = mat2grid(var, x = lon, y = lat)
A GMTgrid object with 1 layers of type Float32
Gridline node registration used
x_min: 0.0      x_max :359.0    x_inc :1.0      n_columns :360
y_min: -89.5    y_max :89.5     y_inc :1.0      n_rows :180
z_min: 0.0      z_max :0.0
Mem layout:     BCB

The horizontal stray lines are due to an unfortunate GMT bug that has proven very elusive. It often disappear when we change the coastlines resolution.

https://github.com/GenericMappingTools/gmt/issues/8235 https://github.com/GenericMappingTools/gmt/issues/3667

melodyjulia commented 1 month ago

If the longitude coordinates are changed to the following:

lon = 0.5:1:359.5

Why you say that?

I found the horizontal stray lines depend on the longitude coordinates, if the coordinates start from 0, such as [0, 1, ...], the stray lines occur. If the coordinates start from 0.5, such as [0.5, 1.5, ...], the stray lines disappear.

joa-quim commented 1 month ago

As I said, the bug (in GMT, not in GMT.jl) is very illusive. For example, now no stray lines.

lon = -180:180; lat = -90:90; var = zeros(Float32,181,361);
G = mat2grid(var, x = lon, y = lat);
viz(G, proj=(name=:Winkel, center=180), coast=true)

nobug

And nor with

lon = 0:360; lat = -90:90; var = zeros(Float32,181,361);
G = mat2grid(var, x = lon, y = lat);
viz(G, proj=:Winkel, coast=true)
melodyjulia commented 1 month ago

I see, thanks!

joa-quim commented 1 month ago

Closing this because it's not a GMT.jl issue