diff --git a/src/GMT.jl b/src/GMT.jl index 2369a8cf1..43d0845f5 100644 --- a/src/GMT.jl +++ b/src/GMT.jl @@ -132,7 +132,7 @@ export stereonet, stereonet!, add_opt, isgeog, numel, scan_opt, extrema_nan, parse_RIr, close_PS_file, getsize, meshgrid, parse_B, parse_BJR, parse_I, - parse_J, parse_R, ressurectGDAL, CPTaliases, isJLL, POSTMAN, + parse_J, parse_R, ressurectGDAL, CPTaliases, isJLL, POSTMAN, PSname, mgd77magref, magref, diff --git a/src/contourf.jl b/src/contourf.jl index cff76ceb9..5f6fae4e1 100644 --- a/src/contourf.jl +++ b/src/contourf.jl @@ -239,6 +239,6 @@ function auto_contour_interval(range) end # --------------------------------------------------------------------------------------------------- -contourf!(cmd0::String="", arg1=nothing, arg2=nothing; kw...) = contourf(cmd0, arg1, arg2, false, KW(kw)) +contourf!(cmd0::String="", arg1=nothing, arg2=nothing; kw...) = contourf(cmd0, arg1, arg2; first=false, kw...) contourf(arg1, arg2=nothing; kw...) = contourf("", arg1, arg2, true, KW(kw)) contourf!(arg1, arg2=nothing; kw...) = contourf("", arg1, arg2, false, KW(kw)) \ No newline at end of file diff --git a/src/drawing.jl b/src/drawing.jl index dea92a05d..d31ce3806 100644 --- a/src/drawing.jl +++ b/src/drawing.jl @@ -83,7 +83,7 @@ function helper_shapes(x, y, cmd; Vd=0, kw...) ((isa(first, Bool) && !first) || (isa(first, Int) && first == 0)) && (CTRLshapes.first[1] = false) if (CTRLshapes.first[1]) - CTRLshapes.fname[1] = GMT.PSname[1] + CTRLshapes.fname[1] = PSname[1] if ((val = find_in_dict(d, [:units])[1]) !== nothing) _cmd, opt_B = parse_B(d, "", DEF_FIG_AXES[1]) cmd *= " --PROJ_LENGTH_UNIT=p" diff --git a/src/gadm.jl b/src/gadm.jl index 9ea7347dd..632cb8487 100644 --- a/src/gadm.jl +++ b/src/gadm.jl @@ -47,7 +47,7 @@ function gadm(country, subregions...; children::Bool=false, names::Bool=false, c reportlevels && return nlayers function _filterlayer(layer, layer_name, level, value, all::Bool=false) - filtered = Vector{GMT.Gdal.Feature}(undef,0) + filtered = Vector{Gdal.Feature}(undef,0) # They Fck changed this between version 3.6 and 4.5. Before, it was always NAME_0, Name_1, ... key = "NAME_$(level)" if (level == 0) diff --git a/src/imgtiles.jl b/src/imgtiles.jl index 6b63f9734..4ee91c9e4 100644 --- a/src/imgtiles.jl +++ b/src/imgtiles.jl @@ -1006,21 +1006,21 @@ function geocoder(address::String; options=String[])::GDtype # Get the geocoder info for a given address. Adapted from https://www.itopen.it/geocoding-with-gdal/ _ops = isempty(options) ? C_NULL : options # The default is ["SERVICE", "OSM_NOMINATIM"] - hSession = GMT.Gdal.OGRGeocodeCreateSession(_ops) - hLayer = GMT.Gdal.OGRGeocode(hSession, address, C_NULL, [""]) - hFDefn = GMT.Gdal.OGR_L_GetLayerDefn(hLayer) - hFeature = GMT.Gdal.OGR_L_GetNextFeature(hLayer) - count = GMT.Gdal.OGR_FD_GetFieldCount(hFDefn) + hSession = Gdal.OGRGeocodeCreateSession(_ops) + hLayer = Gdal.OGRGeocode(hSession, address, C_NULL, [""]) + hFDefn = Gdal.OGR_L_GetLayerDefn(hLayer) + hFeature = Gdal.OGR_L_GetNextFeature(hLayer) + count = Gdal.OGR_FD_GetFieldCount(hFDefn) (count == 0) && (@warn("No result found for the address $address"); return GMTdataset()) dic = Dict{String,String}() for k = 0:count-1 - hFieldDefn = GMT.Gdal.OGR_FD_GetFieldDefn(hFDefn,k) - ((val = GMT.Gdal.OGR_F_GetFieldAsString(hFeature, k)) != "") && (dic[GMT.Gdal.OGR_Fld_GetNameRef(hFieldDefn)] = val) + hFieldDefn = Gdal.OGR_FD_GetFieldDefn(hFDefn,k) + ((val = Gdal.OGR_F_GetFieldAsString(hFeature, k)) != "") && (dic[Gdal.OGR_Fld_GetNameRef(hFieldDefn)] = val) end BB = parse.(Float64, split(dic["boundingbox"], ",")) - GMT.Gdal.OGRGeocodeFreeResult(hLayer) - GMT.Gdal.OGRGeocodeDestroySession(hSession) + Gdal.OGRGeocodeFreeResult(hLayer) + Gdal.OGRGeocodeDestroySession(hSession) D::GDtype = mat2ds([parse(Float64, dic["lon"]) parse(Float64, dic["lat"])], attrib=dic, proj4=prj4WGS84, geom=wkbPoint) D.ds_bbox = [BB[3], BB[4], BB[1], BB[2]] return D diff --git a/src/imshow.jl b/src/imshow.jl index 84bb12c65..0c956c7e9 100644 --- a/src/imshow.jl +++ b/src/imshow.jl @@ -79,8 +79,8 @@ function imshow(arg1, x::AbstractVector{Float64}=Float64[], y::AbstractVector{Fl isFV3D = (isa(arg1, GMTfv) && size(arg1.verts, 2) > 2) end CTRL.limits[1:4] = ginfo[1:4]; CTRL.limits[7:10] = ginfo[1:4] - call_plot3 = ((isa(arg1, GMTdataset) && arg1.geom == Gdal.wkbPolygonZM) || - (isa(arg1, Vector{<:GMTdataset}) && arg1[1].geom == Gdal.wkbPolygonZM) || + call_plot3 = ((isa(arg1, GMTdataset) && arg1.geom == wkbPolygonZM) || + (isa(arg1, Vector{<:GMTdataset}) && arg1[1].geom == wkbPolygonZM) || isFV3D) ? true : false # Should evolve into a fun that detects the several plot3d cases. !call_plot3 && (call_plot3 = isplot3(kw)) return (call_plot3) ? plot3d(arg1; show=see, kw...) : plot(arg1; show=see, kw...) diff --git a/src/maregrams.jl b/src/maregrams.jl index 893f3fad9..cef053302 100644 --- a/src/maregrams.jl +++ b/src/maregrams.jl @@ -112,7 +112,7 @@ function maregrams(x::Real, y::Real; days=2, starttime::String="", printurl::Boo # Find the closest station to input x,y @assert -180 <= x <= 360 && -90 <= y <= 90 "Coordinates must be between -180 and 360 and -90 and 90." d = read_maregrams() - dists = mapproject(d[:pos], G="$x/$y", o=2) + dists = gmt("mapproject -G$(x)/$(y) -o2", d[:pos]) ind = argmin(dists) maregrams(code=d[:code][ind], days=days, starttime=starttime, printurl=printurl) end diff --git a/src/pcolor.jl b/src/pcolor.jl index 31b710978..b8814912d 100644 --- a/src/pcolor.jl +++ b/src/pcolor.jl @@ -42,7 +42,7 @@ This form takes a grid (or the file name of one) as input an paints it's cell wi --- pcolor(GorD; kwargs...) -If `GorD` is either a GMTgrid or a GMTdataset containing a Pearson correlation matrix obtained with ``GMT.cor()``, +If `GorD` is either a GMTgrid or a GMTdataset containing a Pearson correlation matrix obtained with ``cor()``, the processing recieves a special treatment. In this case, other than the `labels` keyword, user is also interested in seing if the automatic choice of x-annotaions angle is correct. If not, one can force it by setting the `rotx` (ot `slanted`) keywords. @@ -73,7 +73,7 @@ by setting the `rotx` (ot `slanted`) keywords. Display a Pearson's correlation matrix ```julia -pcolor(GMT.cor(rand(4,4)), labels=:y, colorbar=1, show=true) +pcolor(cor(rand(4,4)), labels=:y, colorbar=1, show=true) ``` """ function pcolor(X_::VMr, Y_::VMr, C::Union{Nothing, AbstractMatrix{<:Real}}=nothing; first::Bool=true, kwargs...) @@ -155,7 +155,7 @@ function pcolor(cmd0::String="", arg1=nothing; first=true, kwargs...) (cmd0 != "") && (arg1 = gmtread(cmd0)) (isa(arg1, Matrix)) && (arg1 = mat2grid(Float32.(arg1))) isdataframe(arg1) && (arg1 = df2ds(arg1)) - if ((arg1 == arg1' && arg1[1] == 1 && arg1[end] == 1)) # A corr matrix (computed with GMT.cor()) + if ((arg1 == arg1' && arg1[1] == 1 && arg1[end] == 1)) # A corr matrix (computed with cor()) return pcolor(mat2ds(arg1.z); first=first, kwargs...) elseif (isa(arg1, GMTdataset)) # Arrive here when arg1 was originally a DataFrame return pcolor(arg1; first=first, kwargs...) @@ -177,7 +177,7 @@ function pcolor(D::GMTdataset; first=true, kwargs...) end d[:xticks] = (colnames, ang) d[:yticks] = colnames[end:-1:1] - if (D == D' && D[1] == 1 && D[end] == 1) # A correlation matrix (that computed with GMT.cor()) + if (D == D' && D[1] == 1 && D[end] == 1) # A correlation matrix (that computed with cor()) (is_in_dict(d, CPTaliases) === nothing) && (d[:C] = makecpt(T=(-1,1.0,0.1), C="tomato,azure1,dodgerblue4", Z=true)) for n = 1:size(z,2), m = 1:size(z,1) m < n && (z[m,n] = NaN) # Since the matrix is symmetric, remove the upper triangle diff --git a/src/plot.jl b/src/plot.jl index 29ba38025..dc3a15ee4 100644 --- a/src/plot.jl +++ b/src/plot.jl @@ -472,7 +472,7 @@ end function _bar(arg, first::Bool, d::Dict{Symbol,Any}) do_cat = ((haskey(d, :stack) || haskey(d, :stacked)) && isvector(arg) && length(arg) > 2) ? false : true - is_waterfall = ((val = hlp_desnany_str(d, [:stack :stacked])) !== "" && startswith(val, "water")) + is_waterfall = ((val = hlp_desnany_str(d, [:stack :stacked], false)) !== "" && startswith(val, "water")) if (is_waterfall) isa(arg, Vector) && (arg = reshape(arg, 1, length(arg))) # Waterfall stacks must be matrices (arg[1] != 0) && (arg = hcat(repeat([1.0],size(arg,1)), arg)) # If first el != 0 assume coord is missing diff --git a/src/proj_utils.jl b/src/proj_utils.jl index f8f24d371..1e87d7b05 100644 --- a/src/proj_utils.jl +++ b/src/proj_utils.jl @@ -359,7 +359,8 @@ function buffergeo(line::Matrix{<:Real}; width=0, unit=:m, np=120, flatstart=fal # be able to use GMT's mapproject to find only the points that are at the exact 'width' distance from the # line, but at this moment there seems to be an issue in GMT and distances are shorter. However we can # do some cheap statistics to get read of the more inner points a get an almost perfec buffer. - Ddist2line = mapproject(D, L=(line=line, unit=:e)) # Find the distance from buffer points to input line + #Ddist2line = mapproject(D, L=(line=line, unit=:e)) # Find the distance from buffer points to input line + Ddist2line = gmt("mapproject -L+ue", D, line) # Find the distance from buffer points to input line if (isa(Ddist2line, GMTdataset)) d_mean = mean(view(Ddist2line.data, :, 3)) ind = view(Ddist2line.data, :, 3) .>= d_mean @@ -486,7 +487,7 @@ function geodesic_long(lonlat1::VMr, lonlat2::VMr; step=0.0, np=180, proj::Strin dest2, = geod(lonlat1, azA, 40050, unit=:km, proj=proj, epsg=epsg) # Destination of a point after the perimeter t = geodesic([dest1[1] dest1[2]; dest2[1] dest2[2]], step=1000, proj=proj, epsg=epsg) # Temp geodesic arround the perimeter point gmtwrite(_name, t); # Workaround a bug in <= 6.4.0 - x = mapproject([lonlat1[1] lonlat1[2]], L=_name) # Assume the closest point to A corresponds to full perimeter. + x = gmt("mapproject -L$_name", [lonlat1[1] lonlat1[2]]) # Assume the closest point to A corresponds to full perimeter. d, = invgeod(dest1, [x[4] x[5]], proj=proj, epsg=epsg) # Distance from dest1 to the point where the geodesic does a full perimeter. dtot = 39950000 + d # Total length of this geodesic (A bit in excess probably due to geodesics not closing) (step != 0) && (np = round(Int, dtot / step) + 1) # Takes precedence over 'np' because that one has a default value. @@ -494,7 +495,7 @@ function geodesic_long(lonlat1::VMr, lonlat2::VMr; step=0.0, np=180, proj::Strin function get_mindist(az) D = geod(lonlat1, az, linspace(0,dtot,360), proj=proj, epsg=epsg)[1] gmtwrite(_name, D); - mapproject([lonlat2[1] lonlat2[2]], L=_name) # Find closest point to B in the passing by geodesic + gmt("mapproject -L$_name", [lonlat2[1] lonlat2[2]]) # Find closest point to B in the passing by geodesic. end # Using the azimuth computed from the shortest distance between A & B will always give only an aproximation @@ -512,7 +513,7 @@ function geodesic_long(lonlat1::VMr, lonlat2::VMr; step=0.0, np=180, proj::Strin az = (dt < d0) ? az + (d0 / dd) * daz : az - (d0 / dd) * daz d = geod(lonlat1, az, dtot-distAB)[1] # First estimate of the arrival point - D = mapproject(lonlat2, G=(d[1], d[2])) # Distance between estimation and true point + D = gmt("mapproject -G$(d[1])/$(d[2])", (isa(lonlat2, Vector{<:Real})) ? reshape(lonlat2, 1, :) : lonlat2) # Distance between estimation and true point total_len = dtot - distAB - D[end] # Correct 'dtot' that was estimated in excess. p,a = geod(lonlat1, az, linspace(0,total_len,np)) # Possibly because the path is slightly hellicoidal. prj = (proj != "" && (contains(proj,"=lon") || contains(proj,"=lat"))) ? proj : prj4WGS84 diff --git a/src/psxy.jl b/src/psxy.jl index f6b63fccf..c6d36d47b 100644 --- a/src/psxy.jl +++ b/src/psxy.jl @@ -36,7 +36,7 @@ function _common_plot_xyz(cmd0::String, arg1, caller::String, O::Bool, K::Bool, cmd, is_gridtri, arg1 = parse_grid2tri_case(d, cmd, caller, is3D, isFV, O, arg1) isa(arg1, GMTdataset) && (arg1 = with_xyvar(d, arg1)) # See if we have a column request based on column names - if ((val = hlp_desnany_int(d, [:decimate])) !== -999) # Asked for a clever data decimation? + if ((val = hlp_desnany_int(d, [:decimate], del=true)) !== -999) # Asked for a clever data decimation? dec_factor::Int = (val == 1) ? 10 : val arg1 = lttb(arg1, dec_factor) end diff --git a/src/spatial_funs.jl b/src/spatial_funs.jl index 4f9b4165f..bf8e707ff 100644 --- a/src/spatial_funs.jl +++ b/src/spatial_funs.jl @@ -646,7 +646,7 @@ function Base.:in(D1::GMTdataset, D2::GMTdataset)::Union{Bool, Int, Vector{Int}} (!is2Pol && indPoint == 0 && indPol == 0) && error("Input arguments must have a Point and Polygon geometries, or both be Polygons.") if (indPoint == 1 && indPol == 2) inwhichpolygon([D2], D1.data) elseif (indPoint == 2 && indPol == 1) inwhichpolygon([D1], D2.data) - else GMT.contains(D1, D2) # Returns `true` if D1 contains D2. + else contains(D1, D2) # Returns `true` if D1 contains D2. end end diff --git a/src/utils_types.jl b/src/utils_types.jl index 56649bf53..b45187abf 100644 --- a/src/utils_types.jl +++ b/src/utils_types.jl @@ -361,8 +361,8 @@ function _mat2ds(mat::Array{T,N}, txt::Union{String,Vector{String}}, hdr::Vector # By default convert to Doubles, except if instructed to NOT to do it. #(find_in_dict(d, [:datatype])[1] === nothing) && (eltype(mat) != Float64) && (mat = Float64.(mat)) - _geom::Int = Int((geom == 0 && (2 <= length(mat) <= 3)) ? Gdal.wkbPoint : (geom == 0 ? Gdal.wkbUnknown : geom)) # Guess geom - (multi && _geom == 0 && size(mat,1) == 1) && (_geom = Int(Gdal.wkbPoint)) # One row with many columns and MULTI => Points + _geom::Int = Int((geom == 0 && (2 <= length(mat) <= 3)) ? wkbPoint : (geom == 0 ? wkbUnknown : geom)) # Guess geom + (multi && _geom == 0 && size(mat,1) == 1) && (_geom = Int(wkbPoint)) # One row with many columns and MULTI => Points if (isempty(xx)) # No coordinates transmitted if (ndims(mat) == 3) coln = fill_colnames(coln, size(mat,2)-2, is_geog) @@ -959,7 +959,7 @@ function tabletypes2ds(arg, interp=0) #(arg === nothing || isa(arg, GDtype) || isa(arg, Matrix{<:Real})) && return arg isdataframe(arg) && return df2ds(arg) # DataFrames are(?) easier to deal with. isODE(arg) && return ODE2ds(arg, interp=interp) # DifferentialEquations type is a complex beast. - (isa(arg, GMT.Gdal.AbstractDataset) || isa(arg, GMT.Gdal.AbstractGeometry)) && return gd2gmt(arg) + (isa(arg, Gdal.AbstractDataset) || isa(arg, Gdal.AbstractGeometry)) && return gd2gmt(arg) # Guesswork, it may easily screw. colnames = [i for i in fields(arg) if Base.nonmissingtype(eltype(getproperty(arg, i))) <: AbstractFloat]