Add a the leepacifc projection to map the Pacific Ocean with low deformations.

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,

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)

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...)

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
 

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])

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])
 

src/utils_project.jl

@@ -75,49 +75,61 @@ 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")
 
 	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))"))