From 6d248f1973fa048aa89038cadc732b7c1f9a3cf3 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 20 Aug 2025 20:16:36 +0100 Subject: [PATCH 1/5] Add help_desanys --- src/common_options.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/common_options.jl b/src/common_options.jl index 3600fb3d5..9511d0182 100644 --- a/src/common_options.jl +++ b/src/common_options.jl @@ -858,9 +858,9 @@ function parse_proj(p::NamedTuple)::Tuple{String, Bool} # may be absent, but not BOTH, an create a GMT -J syntax string (note: for some projections 'center' # maybe a scalar but the validity of that is not checked here). d = nt2dict(p) # Convert to Dict - if ((val = find_in_dict(d, [:name])[1]) !== nothing) - prj::String, mnemo = parse_proj(string(val)::String) - if (prj != "Cyl_stere" && prj == string(val)::String) + if ((val = hlp_desnany_str(d, [:name])) != "") + prj::String, mnemo = parse_proj(val) + if (prj != "Cyl_stere" && prj == val) @warn("Very likely the projection name ($prj) is unknown to me. Expect troubles") end else @@ -892,7 +892,7 @@ function parse_proj(p::NamedTuple)::Tuple{String, Bool} end # Piggy-back `center`. Things can be wrong here is user does stupid things like using zone & parallels - ((val = find_in_dict(d, [:zone])[1]) !== nothing) && (center = string(val)) + ((val = hlp_desnany_str(d, [:zone])) !== "") && (center = val) if (center == "" && parallels != "") center = "0/0" * parallels elseif (center != "") center *= parallels # even if par == "" @@ -1626,7 +1626,7 @@ function parse_r(d::Dict, cmd::String, del::Bool=true) # Accept both numeric (0 or != 0) and string/symbol arguments opt_val::String = "" if ((val = find_in_dict(d, [:r :reg :registration], del)[1]) !== nothing) - (isa(val, String) || isa(val, Symbol)) && (opt_val = string(" -r",val)[1:4]) + (isa(val, StrSymb)) && (opt_val = string(" -r",val)[1:4]) (isa(val, Integer)) && (opt_val = (val == 0) ? " -rg" : " -rp") cmd *= opt_val end @@ -1656,8 +1656,8 @@ end # --------------------------------------------------------------------------------------------------- function parse_append(d::Dict, cmd::String)::String - if ((val = find_in_dict(d, [:append], true)[1]) !== nothing) - cmd *= " >> " * val::String + if ((val = hlp_desnany_str(d, [:append])) !== "") + cmd *= " >> " * val end return cmd end @@ -1943,7 +1943,7 @@ function opt_pen(d::Dict, opt::Char, symbs::VMs)::String out = string(" -", opt, pen) else if ((val = find_in_dict(d, symbs)[1]) !== nothing) - if (isa(val, String) || isa(val, Real) || isa(val, Symbol)) + if (isa(val, StrSymb) || isa(val, Real)) out = string(" -", opt, val) elseif (isa(val, Tuple)) # Like this it can hold the pen, not extended atts out = string(" -", opt, parse_pen(val)) @@ -2263,7 +2263,7 @@ function prepare2geotif(d::Dict, cmd::Vector{String}, opt_T::String, O::Bool)::T return cmd, opt_T end cmd[1] = helper2geotif(cmd[1]) - if (isa(val, String) || isa(val, Symbol)) # A transparent KML + if (isa(val, StrSymb)) # A transparent KML if (startswith(string(val)::String, "trans")) opt_T = " -TG -W+k" else opt_T = string(" -TG -W+k", val) # Whatever 'val' is end @@ -2287,7 +2287,7 @@ function add_opt_1char(cmd::String, d::Dict, symbs::Vector{Matrix{Symbol}}; del: for opt in symbs ((val = find_in_dict(d, opt, del)[1]) === nothing) && continue args::String = "" - if (isa(val, String) || isa(val, Symbol)) + if (isa(val, StrSymb)) ((args = arg2str(val)) != "") && (args = string(args[1])) elseif (isa(val, Tuple)) for k = 1:numel(val) From db38057471cc2e43a6bb9bcd79d90a6bf91532dd Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 20 Aug 2025 20:17:09 +0100 Subject: [PATCH 2/5] Change var name. --- src/gdal_tools.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gdal_tools.jl b/src/gdal_tools.jl index 1405e1b71..1bcfa03ac 100644 --- a/src/gdal_tools.jl +++ b/src/gdal_tools.jl @@ -360,8 +360,8 @@ function helper_run_GDAL_fun(f::Function, indata, dest::String, opts, method::St if (f == gdalfillnodata!) # This guy has its own peculiarities f(dataset; d...) - o = gd2gmt(dataset) - indata.z, indata.hasnans, indata.layout = o.z, 1, o.layout + D = gd2gmt(dataset) + indata.z, indata.hasnans, indata.layout = D.z, 1, D.layout return nothing elseif (f == gdalsievefilter) o = f(dataset; d...) From c5ae12ef6bf3de1d1f3f03495406e0848178090a Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 20 Aug 2025 20:17:54 +0100 Subject: [PATCH 3/5] Robustify Base.:cat(D1::Vector{<:GMTdataset}, D2::Vector{<:GMTdataset}) --- src/grd_operations.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/grd_operations.jl b/src/grd_operations.jl index c95200a3b..5d17e0e7c 100644 --- a/src/grd_operations.jl +++ b/src/grd_operations.jl @@ -362,6 +362,8 @@ Base.:cat(D1::GMTdataset, D2::GMTdataset) = Base.cat([D1], [D2]) function Base.:cat(D1::Vector{<:GMTdataset}, D2::Vector{<:GMTdataset}) # Concat 2 Vector{GMTdataset}. The important point is that the final 'ds_bbox' field gets set correctly # because plot() uses it to set automatic limits. + isempty(D1) && return D2 + isempty(D2) && return D1 D = vcat(D1, D2) for k = 1:2:length(D1[1].ds_bbox) # Udate the cat'ed ds_bbox. This is crutial D[1].ds_bbox[k] = min(D1[1].ds_bbox[k], D2[1].ds_bbox[k]) From 1a3fb70c6c03ca967d1b24fb5a36845fdf5f3880 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 20 Aug 2025 20:28:49 +0100 Subject: [PATCH 4/5] Save the WKT in meta when available --- src/gdal_utils.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/gdal_utils.jl b/src/gdal_utils.jl index dfa06fd37..39801a6df 100644 --- a/src/gdal_utils.jl +++ b/src/gdal_utils.jl @@ -102,6 +102,7 @@ function gd2gmt(_dataset; band::Int=0, bands=Vector{Int}(), sds::Int=0, pad::Int end hdr = [x_min, x_max, y_min, y_max, z_min, z_max, Float64(!is_grid), x_inc, y_inc] prj = getproj(dataset) + prjwkt = startswith(prj, "PROJCS") ? prj : "" (prj != "" && !startswith(prj, "+proj")) && (prj = toPROJ4(importWKT(prj))) (prj == "") && (prj = seek_wkt_in_gdalinfo(gdalinfo(dataset))) is_tp = (layout == "") # if == "" array is rowmajor and hence transposed @@ -117,7 +118,7 @@ function gd2gmt(_dataset; band::Int=0, bands=Vector{Int}(), sds::Int=0, pad::Int end if (is_grid) #(eltype(mat) == Float64) && (mat = Float32.(mat)) - O = mat2grid(mat; hdr=hdr, v=vvalues, v_unit=vname, proj4=prj, names=desc, is_transposed=is_tp) + O = mat2grid(mat; hdr=hdr, v=vvalues, v_unit=vname, proj4=prj, wkt=prjwkt, names=desc, is_transposed=is_tp) O.layout = (layout == "") ? "TRB" : layout isa(mat, Matrix{<:Complex}) && (append!(O.range, [z_im_min, z_im_max])) # Stick the imaginary part limits in the range if (orig_is_UInt16) @@ -130,7 +131,7 @@ function gd2gmt(_dataset; band::Int=0, bands=Vector{Int}(), sds::Int=0, pad::Int end end else - O = mat2img(mat; hdr=hdr, proj4=prj, noconv=true, names=desc, is_transposed=is_tp) + O = mat2img(mat; hdr=hdr, proj4=prj, wkt=prjwkt, noconv=true, names=desc, is_transposed=is_tp) got_fill_val && (O.nodata = fill_val) O.layout = (layout == "") ? (size(O, 3) == 4 ? "TRBA" : "TRBa") : (layout[end] != 'a' ? layout * "a" : layout) if (n_colors > 0) @@ -1048,6 +1049,7 @@ function xy2lonlat(xy::Matrix{<:Real}, s_srs_=""; s_srs="", t_srs=prj4WGS84) isa(t_srs, Int) && (t_srs = epsg2wkt(t_srs)) (s_srs == "") && error("Must specify at least the source referencing system.") D = ogr2ogr(xy, ["-s_srs", s_srs, "-t_srs", t_srs, "-overwrite"]) + isa(D, Vector) && contains(s_srs, "=lee_os") && length(D) > 1 && return [D[1][1,1] D[1][1,2]; D[2][2,1] D[2][2,2]] # GDAL BUG? return D.data # Return only the array because that's what was sent in end From a7bab28048682207c87fd38eaf68e05eb06c103b Mon Sep 17 00:00:00 2001 From: Joaquim Date: Wed, 20 Aug 2025 20:29:49 +0100 Subject: [PATCH 5/5] Add a the `leepacifc` projection to map the Pacific Ocean with low deformations. --- src/GMT.jl | 2 +- src/grdcut.jl | 1 + src/utils_project.jl | 128 +++++++++++++++++++++++++++++++------------ 3 files changed, 95 insertions(+), 36 deletions(-) diff --git a/src/GMT.jl b/src/GMT.jl index 3ed56f5e8..fad675c2e 100644 --- a/src/GMT.jl +++ b/src/GMT.jl @@ -140,7 +140,7 @@ export blendimg!, lonlat2xy, xy2lonlat, df2ds, mat2ds, mat2grid, mat2img, slicecube, cubeslice, linspace, logspace, fileparts, fields, flipud, fliplr, flipdim, flipdim!, grdinterpolate, pow, tic, toc, theme, tern2cart, geodetic2enu, cpt4dcw, getregion, gd2gmt, gmt2gd, gdalread, gdalshade, gdalwrite, gadm, xyzw2cube, coastlinesproj, graticules, orbits, orbits!, - plotgrid!, worldrectangular, worldrectgrid, + plotgrid!, leepacific, worldrectangular, worldrectgrid, togglemask, earthregions, gridit, grid2tri, magic, rescale, stackgrids, delrows, setgrdminmax!, meshgrid, cart2pol, pol2cart, diff --git a/src/grdcut.jl b/src/grdcut.jl index 6fac3f4ae..5e88ba594 100644 --- a/src/grdcut.jl +++ b/src/grdcut.jl @@ -177,6 +177,7 @@ function axes2pix(xy, dims, x, y, reg=0, layout::String="TC") one_or_zero = (reg == 0) ? 1.0 : 0.0 slope = (dims[col_dim] - one_or_zero) / (x[end] - x[1]); isnan(slope) && (slope = 1.0) # Vertical slices of a cube pix_x = round.(Int, slope .* (xy[:,1] .- x[1]) .+ [1.0, one_or_zero]) + pix_x[1] < 1 && (pix_x[1] = 1) # Happened in worldrectangular() when the grid x limits returned by gdalwarp were not suitable slope = (dims[row_dim] - one_or_zero) / (y[end] - y[1]); isnan(slope) && (slope = 1.0) pix_y = round.(Int, slope .* (xy[:,2] .- y[1]) .+ [1.0, one_or_zero]) diff --git a/src/utils_project.jl b/src/utils_project.jl index 49ca0f74c..8194ae94f 100644 --- a/src/utils_project.jl +++ b/src/utils_project.jl @@ -75,41 +75,53 @@ function worldrectangular(GI::GItype; proj::StrSymb="+proj=vandg +over", pm=0, l _proj = isa(proj, Symbol) ? string(proj) : proj (latlim === nothing && latlims !== nothing) && (latlim = latlims) # To accept both latlim & latlims autolat = false - isa(latlim, StrSymb) && (autolat = true; latlim = nothing) - _latlim = (latlim === nothing) ? (-90,90) : (isa(latlim, Real) ? (-latlim, latlim) : extrema(latlim)) + is_lee_os = contains(_proj, "lee_os") + if (!is_lee_os) + isa(latlim, StrSymb) && (autolat = true; latlim = nothing) + _latlim = (latlim === nothing) ? (-90,90) : (isa(latlim, Real) ? (-latlim, latlim) : extrema(latlim)) + else + _latlim = latlims + end (_latlim[1] == _latlim[2]) && error("Your two triming latitudes ('latlim') are equal.") (_latlim[1] < -90 || _latlim[2] > 90) && error("Don't kid, latlim does not have real latitude limits.") !startswith(_proj, "+proj=") && (_proj = "+proj=" * _proj) - !contains(_proj, " +over") && (_proj *= " +over") + !is_lee_os && !contains(_proj, " +over") && (_proj *= " +over") # lee_os is an exception in this function (pm > 180) && (pm -= 360); (pm < -180) && (pm += 360) - # The pad is an attempt to minimize the empties that may occur in the corners for some projs - pad = (contains(proj, "wintri") && pad == 0) ? 120 : pad # Winkel Tripel needs more - (pad == 0) && (pad = 90) # The default value - res = (isa(coast, Symbol) || isa(coast, String)) ? string(coast) : (coast == 1 ? "crude" : "none") + res = (isa(coast, Symbol) || isa(coast, String)) ? string(coast) : (coast == 1 ? (is_lee_os ? "low" : "crude") : "none") coastlines = (res == "none" && isa(coast, GDtype)) ? coast : GMTdataset{Float64,2}[] # See if we have an costlines argin. (isempty(coastlines) && !isa(coast, Bool) && res[1] != 'c' && res[1] != 'l' && res[1] != 'i' && res[1] != 'h') && error("Bad input for the 'coast' option.") - if (pm >= 0) - Gr = grdcut(GI, R=(-180,-180+pad+pm, -90,90)) # Cut on West to be added at East - Gr.range[1], Gr.range[2] = 180, 180+pad+pm # pretend it is beyound 180 - if (pm < pad) - Gl = grdcut(GI, R=(180-(pad-pm), 180, -90,90)) # Cut from East - Gl.range[1], Gl.range[2] = Gl.range[1]-360., Gl.range[2]-360. # To be added on the West + if (!is_lee_os) + # The pad is an attempt to minimize the empties that may occur in the corners for some projs + pad = (contains(proj, "wintri") && pad == 0) ? 120 : pad # Winkel Tripel needs more + (pad == 0) && (pad = 90) # The default value + if (pm >= 0) + Gr = grdcut(GI, R=(-180,-180+pad+pm, -90,90)) # Cut on West to be added at East + Gr.range[1], Gr.range[2] = 180, 180+pad+pm # pretend it is beyound 180 + if (pm < pad) + Gl = grdcut(GI, R=(180-(pad-pm), 180, -90,90)) # Cut from East + Gl.range[1], Gl.range[2] = Gl.range[1]-360., Gl.range[2]-360. # To be added on the West + else + Gl = grdcut(GI, R=(-180+pm-pad, 180, -90,90)) # Trim original at the West + end else - Gl = grdcut(GI, R=(-180+pm-pad, 180, -90,90)) # Trim original at the West + Gl = grdcut(GI, R=(180-(pad-pm), 180, -90,90)) # Cut on East to be added at West + Gl.range[1], Gl.range[2] = -180-(pad-pm), -180 # pretend it is beyound -180 + if (pm > -pad) + Gr = grdcut(GI, R=(-180, -180+(pad+pm), -90,90)) # Cut from West + Gr.range[1], Gr.range[2] = Gr.range[1]+360., Gr.range[2]+360. # To be added on the East + else + Gr = grdcut(GI, R=(-180, 180+pad+pm, -90,90)) # Trim original at the East + end end - else - Gl = grdcut(GI, R=(180-(pad-pm), 180, -90,90)) # Cut on East to be added at West - Gl.range[1], Gl.range[2] = -180-(pad-pm), -180 # pretend it is beyound -180 - if (pm > -pad) - Gr = grdcut(GI, R=(-180, -180+(pad+pm), -90,90)) # Cut from West - Gr.range[1], Gr.range[2] = Gr.range[1]+360., Gr.range[2]+360. # To be added on the East - else - Gr = grdcut(GI, R=(-180, 180+pad+pm, -90,90)) # Trim original at the East + if (pm >= 0 && pm < pad) || (pm < 0 && pm > -pad) G = grdpaste(Gl,GI); G = grdpaste(G,Gr) + else G = grdpaste(Gl,Gr) end - end - if (pm >= 0 && pm < pad) || (pm < 0 && pm > -pad) G = grdpaste(Gl,GI); G = grdpaste(G,Gr) - else G = grdpaste(Gl,Gr) + west, east = -180.0, 180.0 + else + G = GI + #west, east = 90.0, 300.0 + west, east = G.range[1]+20, G.range[2]-30 # Use the range of the input grid, that should ... end (pm != 0) && (_proj *= " +pm=$pm") @@ -117,7 +129,7 @@ function worldrectangular(GI::GItype; proj::StrSymb="+proj=vandg +over", pm=0, l lims_geog = G.range G = gdalwarp(G, ["-t_srs", _proj]) - xy = lonlat2xy([-180.0+pm 0; 180+pm 0], t_srs=_proj) + xy = lonlat2xy([west+pm 0; east+pm 0], t_srs=_proj) pix_x = axes2pix(xy, size(G), [G.x[1], G.x[end]], [G.y[1], G.y[end]], G.registration, G.layout)[1] if (autolat) t = G[pix_x[1], :] # The column corresponding to lon = -180 @@ -133,12 +145,56 @@ function worldrectangular(GI::GItype; proj::StrSymb="+proj=vandg +over", pm=0, l G = grdcut(G, R=(G.x[pix_x[1]], G.x[pix_x[2]], yb, yt)) G.remark = "pad=$pad" # Use this as a pocket to use in worldrectgrid() - return (res != "none" || !isempty(coastlines)) ? (G, worldrectcoast(proj=_proj, res=res, coastlines=coastlines, limits=lims_geog)) : G + if (is_lee_os) + cl = (isempty(coastlines)) ? pscoast(dump=:true, res=res, region=lims_geog) : coastlines + cl = lonlat2xy(cl, t_srs=_proj) + else + cl = worldrectcoast(proj=_proj, res=res, coastlines=coastlines, limits=lims_geog) + end + return (res != "none" || !isempty(coastlines)) ? (G, cl) : G +end + +# ----------------------------------------------------------------------------------------------- +""" + GI[,coast] = leepacific(fname::String; latlims=nothing, lonlims=nothing, coast=true) + +Project a geographical grid/image in the Lee Oblated Stereographic projection centered on the Pacific Ocean. + +- `GI`: A GMTgrid or GMTimage data type. `GI` can also be a string with a file name of a grid or image. + NOTE: This grid/image should have longitudes covering the range 90 to 300 degrees (here lon in [0 360] is better) +- `latlims`: Latitudes used in `region` when reading the grid/image. The default is (-90, 75). +- `lonlims`: Longitudes used in `region` when reading the grid/image. You cannot deviate mutch from (90,300). +- `coast`: Return also the coastlines projected with `proj`. Pass `coast=res`, where `res` is one of + GMT coastline resolutions (*e.g.* :crude, :low, :intermediate). `coast=true` is <==> `coast=:crude` + Pass `coast=D`, where `D` is vector of GMTdataset containing coastline polygons with a provenience + other than the GSHHG GMT database. If `coast=false` the funtion returns only the projected grid/image. + + +### Returns +A grid or an image and optionally the coastlines ... or errors. Not many projections support the procedure +implemented in this function. +The working or not is controlled by PROJ's `+over` option https://proj.org/usage/projections.html#longitude-wrapping + +### Example: + G,cl = leepacific("@earth_relief_10m_g") + grdimage(G, shade=true, plot=(data=cl,), cmap=:geo, B=:none) + plotgrid!(G, grid, show=true) +""" +function leepacific(fname::String; region=(90.0, 300.0, -90.0, 75.0), latlims=nothing, lonlims=nothing, coast=true) + reg = Float64.(collect(region)) + lonlims !== nothing && (reg[1] = lonlims[1]; reg[2] = lonlims[2]) # Use the lonlims if given + latlims !== nothing && (reg[3] = latlims[1]; reg[4] = latlims[2]) # Use the lonlims if given + while(reg[1] < 0) reg[1] += 180.0; reg[2] += 180.0 end + G = gmtread(fname, R=[reg[1]-20.0, reg[2]+30.0, reg[3], reg[4]]) + worldrectangular(G, proj="+proj=lee_os", latlims=(latlims===nothing) ? (-90,75) : latlims, coast=coast) +end +function leepacific(GI::GItype; latlims=nothing, coast=true) + worldrectangular(GI, proj="+proj=lee_os", latlims=(latlims===nothing) ? (-90,75) : latlims, coast=coast) end # ----------------------------------------------------------------------------------------------- """ - cl = coastlinesproj(proj="?", res="crude", coastlines=nothing) + cl = coastlinesproj(proj="?", res="crude", coastlines=nothing, limits=Float64[]) Extract the coastlines from GMT's GSHHG database and project them using PROJ (NOT the GMT projection machinery). This allows the use of many of the PROJ projections that are not available from pure GMT. @@ -151,7 +207,7 @@ This allows the use of many of the PROJ projections that are not available from coastlines that expand more than 360 degrees (`worldrectangular` uses this). ### Returns -A Vector of GMTdataset containing the projected (or not) world GSHHG coastlines at resolution `res`. +A Vector of GMTdataset containing the projected world GSHHG coastlines at resolution `res`. ### Example cl = coastlinesproj(proj="+proj=ob_tran +o_proj=moll +o_lon_p=40 +o_lat_p=50 +lon_0=60"); @@ -188,7 +244,6 @@ function worldrectcoast(; proj::StrSymb="", res="crude", coastlines::Vector{<:GM cl = (isempty(coastlines)) ? coast(dump=:true, res=res, region=length(limits)==4 ? limits : :global) : coastlines (_proj == "" && isempty(limits)) && return cl # No proj required nor extending the coastlines (round || isempty(limits)) && return lonlat2xy(cl, t_srs=_proj) # No extensiom so we are donne. - (_proj != "" && !contains(_proj, " +over")) && (_proj *= " +over") pm = (contains(_proj, "+pm=")) ? parse(Float64, string(split(split(_proj, "+pm=")[2])[1])) : 0.0 pad = !isempty(limits) ? ((limits[2] - limits[1]) - 360.0) / 2 : 0.9 @@ -272,15 +327,15 @@ the `GI` metadata. A Vector of GMTdataset containing the projected meridians and parallels. `grat[i]` attributes store information about that element lon,lat. """ -function worldrectgrid(G_I::GItype; width=(30,20), grid=Vector{Vector{Real}}[], annot_x::Union{Vector{Int},Vector{Float64}}=Int[]) +function worldrectgrid(G_I::GItype; width=(30,20), grid=Vector{Vector{Real}}[], annot_x::Union{Vector{Int},Vector{Float64}}=Int[], worldrect=true) prj = getproj(G_I, proj4=true) pad = contains(G_I.remark, "pad=") ? parse(Int,G_I.remark[5:end]) : 60 - worldrectgrid(proj=prj, width=width, pad=pad, grid=grid, annot_x=annot_x) + worldrectgrid(proj=prj, width=width, pad=pad, grid=grid, annot_x=annot_x, worldrect=worldrect) end -function worldrectgrid(D::GDtype; width=(30,20), grid=Vector{Vector{Real}}[], annot_x::Union{Vector{Int},Vector{Float64}}=Int[]) +function worldrectgrid(D::GDtype; width=(30,20), grid=Vector{Vector{Real}}[], annot_x::Union{Vector{Int},Vector{Float64}}=Int[], worldrect=true) # Normally called only by graticules, not directly prj = getproj(D, proj4=true) - worldrectgrid(proj=prj, width=width, grid=grid, annot_x=annot_x) + worldrectgrid(proj=prj, width=width, grid=grid, annot_x=annot_x, worldrect=worldrect) end function worldrectgrid(; proj::StrSymb="", width=(30,20), grid=Vector{Vector{Real}}[], annot_x::Union{Vector{Int},Vector{Float64}}=Int[], pm=0.0, worldrect=true, pad=60) @@ -290,6 +345,8 @@ function worldrectgrid(; proj::StrSymb="", width=(30,20), grid=Vector{Vector{Rea _proj = isa(proj, Symbol) ? string(proj) : proj _proj == "" && error("Input has no projection info") + is_lee_os = contains(_proj, "lee_os") + is_lee_os && (worldrect = false) !startswith(_proj, "+proj=") && (_proj = "+proj=" * _proj) (contains(_proj, "+pm=")) && (pm = parse(Float64, string(split(split(proj, "+pm=")[2])[1]))) (pm != 0 && !contains(_proj, "+pm=")) && (_proj *= " +pm=$pm") @@ -341,8 +398,9 @@ function worldrectgrid(; proj::StrSymb="", width=(30,20), grid=Vector{Vector{Rea end for p = parallels Dgrid[n+=1] = mat2ds(lonlat2xy([parallel fill(p, length(parallel))], t_srs=_proj), attrib=Dict("para_b" => "$p,$(parallel[1])", "para_e" => "$p,$(parallel[end])", "annot" => "n", "n_meridians" => "$(length(meridians))", "n_parallels" => "$(length(parallels))")) + is_lee_os && 0 <= p <= 20 && (Dgrid[n].data[30:end-24,1] .= NaN) # lee_os proj has problems at low lats. Line reenters domain. end - !worldrect && check_gaps(Dgrid, length(meridians)+1, length(Dgrid)) # Try this only with non-worldrect + !is_lee_os && !worldrect && check_gaps(Dgrid, length(meridians)+1, length(Dgrid)) # Try this only with non-worldrect else # Cartesian graticules for m = meridians Dgrid[n+=1] = mat2ds([fill(m,length(meridian)) meridian], attrib=Dict("merid_b" => "$m,-90", "merid_e" => "$m,90", "n_meridians" => "$(length(meridians))", "n_parallels" => "$(length(parallels))"))