Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

some Longitude regions will cause strange lines #398

Open
htyeim opened this issue Jul 16, 2020 · 5 comments
Open

some Longitude regions will cause strange lines #398

htyeim opened this issue Jul 16, 2020 · 5 comments

Comments

@htyeim
Copy link
Contributor

htyeim commented Jul 16, 2020

If I want to plot a region with longitude from 170 to 140 (which means 170-180 0-140, all without 140-170).
I can't give the region like this (170,140,-90,90), as GMT.jl/gmt will rise a error basemap [ERROR]: Option -R parsing failure. . But I found that the region could be (170-360,140,-90,90) or (170,140+360,-90,90). In most case it work fine, but for large regions like (170-360,155,-90,90) there will be some strange lines...

The code that can produce the strange figures:

using GMT
function test(rlon=[170,128])
    gmtbegin()
    gmtfig(string("test_R_", rlon[1], "_", rlon[2]), fmt="png")
    if rlon[1] > rlon[2]
        rlon[1] -= 360
    end
    baseregion = (rlon[1], rlon[2], -90, 90)
    lon = (baseregion[1] + baseregion[2]) / 2.0
    proj = (name = :Robinson,
            # center = [mean(baseregion[1:2]) mean(baseregion[3:4])],
            center = lon,)
    basemap(region=baseregion, proj=proj, frame=(axes = :WSen,),
                figsize=12, )
    topo = makecpt("-H", color=:geo, range=(-8800, 8800,), continuous=true);

    grdimage!("@earth_relief_20m.grd", region=baseregion,
                nan_alpha=true,
                color=topo,
                colorbar=true)
    coast(region=baseregion, 
            borders=(type = :a, pen = (0.06, :black, "--")),
            shore=(type = 1, pen = (0.6, :black,)), 
            # Vd=1,
            )
    gmtend()
end

for lon_right = 5:10:170
    test([170,lon_right]) 
end

(170,5,-90,90)
test_R_170_5

(170,145,-90,90)
test_R_170_145

(170,155,-90,90)
test_R_170_155

(170,165,-90,90)
test_R_170_165

It seems starts from (170,151) the strange lines appear. I don't know if this is related to #357

@joa-quim
Copy link
Member

It's a GMT bug. We have had bugs of this type since immemorial times. Can be reproduced with

coast(proj=(name=:Robinson, center=-18), region=(-190,155,-90,90),  show=1)

@joa-quim
Copy link
Member

joa-quim commented Jul 16, 2020

However the same command will work if we use a higher coastlines resolution.

coast(proj=(name=:Robinson, center=-18), region=(-190,155,-90,90), resolution=:low, show=1)

when resolution is not set the program uses resolution=automatic, which for these limits turn out to be resolution=:crude. However, although this may solve the current case it may very well show up in other limits combination. I will open an issue in GMT (3667).

@htyeim
Copy link
Contributor Author

htyeim commented Jul 16, 2020

I remember that when I used matplotlib (with cartopy) there are also some strange lines, especially when orthographic projection is used. They are both almost pure python packages but when I try to check the source code it's so hard to read and track variables (maybe I don't have enough knowledge...). Anyway I have changed my main program language to Julia and so far so good.

The main problem I think is related to the geographic longitude range. 180/-180 (0/360) are the same line. When I try to plot large scale maps, I always need to check again and again...

Also I think that it would be great if GMT could accept region range like (90, 60, -90,90) (lon1>lon2) and then plot longitude range from 90 to 180 then from -180 to 60.

@joa-quim
Copy link
Member

joa-quim commented Jul 16, 2020

The problem is that when a polygon is cut in, say, the left side and continues in the right side, that polygon has to be split in two and one of the halves needs to modified to plot correctly on the other side. And this is not only restricted to geographical. What has happen is that there are many combinations depending on the projection and some cases continue to fail. Only Paul has the patience (and knowledge) to keep tracing these failures.

Several years ago we have had a discussion about the possibility of letting w > e in limits but, and I don't remember why, we decided against. Note that this would make sense only in geographical coordinates.

Edit: no actually the discussion was about letting w = e (for some reason).
You may open a feature request in GMT. It would not be too difficult to do this at GMT.jl level but there I cannot check if the coordinates are geographic.

@htyeim
Copy link
Contributor Author

htyeim commented Jul 17, 2020

Thank you.... More complex than I though 😸 I think I can just check the region before I call GMT functions...

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants