option DCW issue

[liming-he]

Might not be a bug. For the following code, if I add "DCW=(country=:CA, pen=(0.225,:black))", the Canada border will be shown, but the grid line (lat, lon) will not show.

If I remove "DCW=(country=:CA, pen=(0.225,:black))", the grid line shows well. Any reason why the country and grid line repel each other?

Thanks. Liming.

for yy = 1: 1
file_out = @sprintf("Pixel_%s_%04d.png", this_variable, year(year_list[yy]));

file_out = joinpath(data_path1, file_out);
println(file_out)

x.z[:] = reverse(transpose(data[:,:,yy]), dims=1);
x.z[x.z .== 0] .= NaN


# figure 7
topo = makecpt(color=:jet, range= (0, 250, 10), continuous=true) #make colorbar  other options: color=:polar
gmt("set MAP_GRID_PEN_PRIMARY thinnest,.");

grdimage(x, show=false, 
    interp = "n", 
    frame=(axes=(:left_bare,:bottom_bare,:right_bare,:top_bare),  annot=true, ),
    coast=(**DCW=(country=:CA, pen=(0.225,:black)),** frame=(axes="WSEN", grid=10,  annot=10, ticks=2, ), ),
    colorbar=(pos=(anchor=:BC,length=(10.5,0.5), horizontal=true, offset=(0,.8), ), color=topo, show=true, frame=(ylabel=:"mm year@+-1@+",ticks = :auto, annot = :auto) ), 
    par=(MAP_ANNOT_ORTHO=:n, MAP_TITLE_OFFSET = 0.8, ), fmt=:png,
    savefig = file_out, dpi=600);

end

[joa-quim]

Hmm, can you please provide an example that I can reproduce?

[liming-he]

here is testing code.

my intention is to plot figure with both grid lines and CA border.

a file with projection is attached (tiff but zipped).

thanks a lot.

Mascons_JPL_result.zip

using GMT

testing_path = "/mnt/hdd1/GRACE_data/Mascons_JPL";
data_path1 = "/mnt/ssd1/Forcing/EALCO_Output";

input

file1 = joinpath(testing_path, "Mascons_JPL_result.Tiff");
x = gmtread(file1, grid=true, band=:1);

output

file_out = joinpath(data_path1, "test1.png");
x.z[:] .= 1:960*1140;

############### version 1 ####################

this produces figure with Canada border but without grid lines.

topo = makecpt(color=:jet, range= (0, 960*1140, 10), continuous=true) #make colorbar other options: color=:polar
gmt("set MAP_GRID_PEN_PRIMARY thinnest,.");
grdimage(x, show=false,
interp = "n",
frame=(axes=(:left_bare,:bottom_bare,:right_bare,:top_bare), annot=true, ),
coast=(DCW=(country=:CA, pen=(0.225,:black)), frame=(axes="WSEN", grid=10, annot=10, ticks=2, ), ),
colorbar=(pos=(anchor=:BC,length=(10.5,0.5), horizontal=true, offset=(0,.8), ), color=topo, show=true, frame=(ylabel=:"GRID ID",ticks = :auto, annot = :auto) ),
par=(MAP_ANNOT_ORTHO=:n, MAP_TITLE_OFFSET = 0.8, ), fmt=:png,
savefig = file_out, dpi=600);

############# version 2 ###############################

this produces figure with shorelines and grid lines.

file_out = joinpath(data_path1, "test2.png");
gmt("set MAP_GRID_PEN_PRIMARY thinnest,.");
topo = makecpt(color=:jet, range= (0, 960*1140, 10), continuous=true) #make colorbar other options: color=:polar
grdimage(x, show=false,
interp = "n",
frame=(axes=(:left_bare,:bottom_bare,:right_bare,:top_bare), annot=true, ),
coast=(frame=(axes="WSEN", grid=10, annot=10, ticks=2, ), ),
colorbar=(pos=(anchor=:BC,length=(10.5,0.5), horizontal=true, offset=(0,.8), ), color=topo, show=true, frame=(ylabel=:"GRID ID",ticks = :auto, annot = :auto) ),
par=(MAP_ANNOT_ORTHO=:n, MAP_TITLE_OFFSET = 0.8, ), fmt=:png,
savefig = file_out, dpi=600);

[liming-he]

[joa-quim]

Sorry, there seem to be issues here but I cant even reproduce your figures. This is what I get

(which GMT.jl version are you using?)

grdimage(x,  interp="n", frame=(axes=(:left_bare,:bottom_bare,:right_bare,:top_bare), annot=true), coast=(frame=(axes="WSEN", grid=10, annot=10, ticks=2)), colorbar=(pos=(anchor=:BC,length=(10.5,0.5), horizontal=true, offset=(0,.8)), color=topo, show=true, frame=(ylabel=:"GRID ID",ticks = :auto, annot = :auto)), par=(MAP_ANNOT_ORTHO=:n, MAP_TITLE_OFFSET=0.8), show=1)

[liming-he]

Yours read the image file correctly. It is strange to see a different version of grid lines.

Mine:
GMT v1.6
Julia 1.9
Ubuntu 20.04

That single line produced below:

[joa-quim]

I don't have an answer for now. This is a harder case where one have a projected grid and want to overlay coastlines that guess the right projection. This is a case that I've worked a couple times (fixing previous failures) but current version (1.10) seems to be behaving worst than the one you have for this case. I will need more time to understand what/why things are failing the way they are.

[liming-he]

Thanks a lot!

Also to mention, I encountered these warning/errors when using grdimage. But the plot looks nice anyway regardless of these errors.

┌ Warning: Que exagero de cores
└ @ GMT ~/.julia/packages/GMT/Mxc9g/src/gmt_main.jl:1082
ERROR 1: Point outside of projection domain
ERROR 10: Pointer 'hTransform' is NULL in 'OCTTransform'.

ERROR 10: Pointer 'hTransform' is NULL in 'OCTTransform'.

[joa-quim]

First warning (in portuguese, sorry) is due to the fact that your makecpt command generates range= (0, 960*1140, 10) 109440 colors. The others are more complex to figure out but are coming from GDAL complaining that data points are outside domain (part of the problem that needs more understanding).

[joa-quim]

OK, I found the regression and now I'm also able to reproduce your results. The issue seems ti lie deep in the GMT (C lib) interaction with GDAL. If you have a standalone GMT CLI, you'll see that this works (it uses the GMT syntax for establishing the projection)

gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -JL-95/0/49/77/15c -Bpa10f2g10 -BWSEN -Da -png lixo

but this one, that uses the PROJ4 syntax, prints many complains and ends up not plotting the grid lines.

gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -J+proj=lcc+lat_0=0+lon_0=-95+lat_1=49+lat_2=77+width=15c/0 -Bpa10f2g10 -BWSEN -Da -png lixo
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
...

[liming-he]

Thanks. Could you share which version of GMT CLI you used?

$ gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -J+proj=lcc+lat_0=0+lon_0=-95+lat_1=49+lat_2=77+width=15c/0 -Bpa10f2g10 -BWSEN -Da -png lixo
gmt [ERROR]: Syntax error: Unrecognized keyword FONT_SUBTITLE.
gmt [ERROR]: Syntax error: GMT_VERBOSE given illegal value (warning)!
gmt [ERROR]: Syntax error: Unrecognized keyword MAP_EMBELLISHMENT_MODE.
gmt [ERROR]: Syntax error: MAP_FRAME_AXES given illegal value (auto)!
gmt [ERROR]: Syntax error: MAP_FRAME_WIDTH given illegal value (autc)!
gmt [ERROR]: Pen attributes not using just - and . for dashes and dots. Offending character --> :

$ gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -JL-95/0/49/77/15c -Bpa10f2g10 -BWSEN -Da -png lixo
gmt [ERROR]: Syntax error: Unrecognized keyword FONT_SUBTITLE.
gmt [ERROR]: Syntax error: GMT_VERBOSE given illegal value (warning)!
gmt [ERROR]: Syntax error: Unrecognized keyword MAP_EMBELLISHMENT_MODE.
gmt [ERROR]: Syntax error: MAP_FRAME_AXES given illegal value (auto)!
gmt [ERROR]: Syntax error: MAP_FRAME_WIDTH given illegal value (autc)!
gmt [ERROR]: Pen attributes not using just - and . for dashes and dots. Offending character --> :

[joa-quim]

Hmm, this makes think you are using an old GMT version

gmt [ERROR]: Syntax error: Unrecognized keyword MAP_EMBELLISHMENT_MODE.

what does this print?

gmt --version

Anyway the Julia wrapper install its own GMT version and that should be GMT6.5.0. The issue you reported roots on GMT itself and I fixed a couple days ago on the master (GMT's, not GMT.jl) version and tried to rebuild the GMT_jl package but the Yggdrasil built is now failing on MacOS and I have no idea why.

[liming-he]

Thanks for helping, - involving to fix problems upstream. I will test it using new version.

name = "GMT"
uuid = "5752ebe1-31b9-557e-87aa-f909b540aa54"
authors = ["Joaquim Luis [email protected]"]
version = "1.6.0"

Binary program via Ubuntu 20.04:
GMT 6.0 (in 2019)

[joa-quim]

No, I'm afraid it's not that simple. You cannot test it with neither of the above. Options are:

1- A new GMT_jll build, but that needs understanding/getting support to know why now BinaryBuild fails on Macs

2- Or build GMT yourself from master and follow instructions in https://github.com/GenericMappingTools/GMT.jl/blob/master/src/GMT.jl#L31

[liming-he]

I am able to see the difference between these two lines with gmt v6.4 in another machine .

gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -JL-95/0/49/77/15c -Bpa10f2g10 -BWSEN -Da -png lixo

gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -J+proj=lcc+lat_0=0+lon_0=-95+lat_1=49+lat_2=77+width=15c/0 -Bpa10f2g10 -BWSEN -Da -png lixo

Before a bug can be fixed, is there a way to pass "-JL-95/0/49/77/15c" into grdimage (coast) in gmt.jl? For example: adding proj="L-95/0/49/77/15c" into coast:

grdimage(x, show=false, 
interp = "n", 
frame=(axes=(:left_bare,:bottom_bare,:right_bare,:top_bare),  annot=true, ),
coast=(proj="L-95/0/49/77/15c", DCW=:CA,  frame=(axes="WSEN", grid=10, annot=10, ticks=2),show=false),
colorbar=(pos=(anchor=:BC,length=(10.5,0.5), horizontal=true, offset=(0,.8), ), color=topo, show=true, frame=(ylabel=Symbol(y_label), ticks = :auto, annot = :auto) ), 
par=(MAP_TITLE_OFFSET = 0.8, ), fmt=:png,
savefig = file_out, dpi=600);

However, this did not fix the problem (or see difference except the errors are gone).

julia> x
title: Grid imported via GDAL
Gridline node registration used
x_min: -2.595e6 x_max :3.095e6 x_inc :10000.0 n_columns :570
y_min: 5.705e6 y_max :1.0495e7 y_inc :10000.0 n_rows :480
z_min: 4.999104976654053 z_max :9999.0
PROJ: +proj=lcc +lat_0=0 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
480×570 Matrix{Float32}:

should I also change the property of x (x.PROJ)? Is there a way to modify the value of x.PROJ in julia, manually?

Thanks.

[joa-quim]

If you update your GMT.jl all of these problems should be solved now.

gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -JL-95/0/49/77/15c -Bpa10f2g10 -BWSEN -Da -png lixo
.
gmt coast -ECA+p0.225,black -R-122.913612/36.214964/-9.977657/62.430681+r -J+proj=lcc+lat_0=0+lon_0=-95+lat_1=49+lat_2=77+width=15c/0 -Bpa10f2g10 -BWSEN -Da -png lixo

These two now give me the same plot with the foxed GMT (C), the same you get when you update GMT.jl

Before a bug can be fixed, is there a way to pass "-JL-95/0/49/77/15c" into grdimage (coast) in gmt.jl? For example: adding proj="L-95/0/49/77/15c" into coast:

Yes but do not add the figure size. That is, this works proj="L-95/0/49/77"

grdimage(x, show=false, 

No need to use show=false. That is the default for all modules except viz|imshow

should I also change the property of x (x.PROJ)? Is there a way to modify the value of x.PROJ in julia, manually?

No, but you can if you want. x is a struct with the fields described in the GMTgrid docs, so you can do
x.proj4 = "+proj=longlat" or just skip the "+proj=" part that should be added automatically when needed (but not 100% sure it happens all the times that is should)

[liming-he]

Thanks.

I used "] add https://github.com/GenericMappingTools/GMT.jl" to add the latest GMT.jl. Now I was able to plot correctly with both grid lines and boundary from DCW (See attached), although some other packages could not be compiled.

Two more questions here:

1. the grid and polygons do not align well for the northern islands. But the grid looks aligned well with polygons in ArcMap.
   

2. the sizes of x and y for a grid are one more than the size(grid.z). This does not look consistent with:
   https://www.generic-mapping-tools.org/GMT.jl/dev/types/

Is that designed so?

file1 = joinpath(testing_path, "can-drainage-region-5km-lcc_tif.tif");
x = gmtread(file1, grid=true);
julia> x
A GMTgrid object with 1 layers of type Float32
title: Grid imported via GDAL
Pixel node registration used
x_min: -2.6e6 x_max :3.1e6 x_inc :10000.0 n_columns :570
y_min: 5.7e6 y_max :1.05e7 y_inc :10000.0 n_rows :480
z_min: 4.999104976654053 z_max :9999.0
Mem layout: BCB
PROJ: +proj=lcc +lat_0=0 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +a=6378137 +rf=298.257024882273 +units=m +no_defs

julia> size(x.x)
(571,)

julia> size(x.y)
(481,)

julia> size(x.z)
(480, 570)

[joa-quim]

I used "] add https://github.com/GenericMappingTools/GMT.jl" to add the latest GMT.jl.

The recommended way is to do ] up GMT or just ] up to update all.

2. the sizes of x and y for a grid are one more than the size(grid.z). This does not look consistent with:

Yes they do and that part is correct. Please see the pixel vs grid registration at https://docs.generic-mapping-tools.org/dev/reference/options.html#grid-registration-the-r-option

1. the grid and polygons do not align well for the northern islands.

Yes, I can confirm that. This type of plots are more complicated than they look at first. Here the grid is already projected and has to be plotted with a linear projection, but to have the frame in geographics one has to compute the right limits (-R) to pass to the basemap module so that it can plot the frame, annotations and graticules correctly. Is this part that is failing to be perfect. Will try to find out why, but meanwhile I think the way to go is to reproject your grid back to geographic (with grdproject or gdalwarp) and do the figure using the Lambert Conformal projection.

[liming-he]

Thank you very much. These are very helpful/useful!

[joa-quim]

What I can tell you is that is the misfit is not due to a projection conversion error. If I reproject the border vector and overlay it with the plot command, I can get a figure like this

My (current) strong suspicion is that the misfit is due to the fact that grid (or the image in my case) central longitude is not equal to the projections lon_0 (in my case that difference is of 0.5 degrees) an that causes a slight rotation of the image that results in a wrong estimation of the -R parameter that is passed (under the hood) to the coast command.

[liming-he]

Any quick fix if this is the case: "that causes a slight rotation of the image that results in a wrong estimation of the -R parameter that is passed (under the hood) to the coast command"? Thanks.

[joa-quim]

That suspicion was not confirmed. I have troubles understanding why (and where) this is failing. The fix should be to reproject the grid to geogs and redo the plot.

[liming-he]

Update:

1. I plot my grid file in ArcMap, it aligns with Canada boundaries well (double confirmed).

2. I plot ("grdimage") the same grid with "+proj=lcc +lat_0=0 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs" imported from "x = gmtread(file1, grid=true);"; that is the first image (misalignment).
   

3. I manually modified x.proj4 = "+proj=lcc +lat_0=0.091 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs" and replot using grdimage, the alignment is okay. Note here the change is from "+lat_0=0" to "+lat_0=0.091".
   

4. with the original x.proj4, I manually added 15000m to array x.y (y axis), and also added 15000m to x.range[3] and x.range[4]; the alignment is okay.
   

Note 1: for 3 and 4, the there is still small misalignment in the very west and east sides.
Note 2: the "15000 m" in y axis and adding 0.091 to lat_0 do not look equal in distance. Maybe my eye could not identify small changes.

[liming-he]

wondering if an example can be provided:
grdimage to plot grid data, user provided shape file, and lat,lon grid lines? or can "coast" be used to provide user's shapefile? thanks.

[joa-quim]

wkt2proj(getproj("Ilcc.tiff"))
"+proj=lcc +lat_0=0 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"

# Reproject to geogs. Note, for some reason (date line crossing?) I had to help it by providing output limits.
Ican_geo = gdalwarp("Ilcc.tiff", ["-t_srs", "+proj=lonlat", "-te", "-138", "40", "-51", "84"]);

viz(Ican_geo, proj="+proj=lcc +lat_0=0 +lon_0=-95 +lat_1=49 +lat_2=77 +x_0=0 +y_0=0 +datum=WGS84", coast=(DCW=(country=:CA, pen=(0.1,:red)),))

[joa-quim]

or can "coast" be used to provide user's shapefile?

The DCW polygons yes. The GSHHG vectors also but they are not contiguous.

D = coast(DCW=:CA, dump=true);
gmtwrite("Canada.shp", D)