GDAL related functions

gdaltranslate(indata, options=String[]; dest="/vsimem/tmp", kwargs...)

Convert raster data between different formats and other operations also provided by the GDAL 'gdal_translate' tool. Namely sub-region extraction and resampling. The kwargs options accept the GMT region (-R), increment (-I), target SRS (-J). Any of the keywords outgrid, outfile or save = outputname options to make this function save the result in disk in the file 'outputname'. The file format is picked from the 'outputname' file extension. When no output file name is provided it returns a GMT object (either a grid or an image, depending on the input type). To force the return of a GDAL dataset use the option gdataset=true.

  • indata: Input data. It can be a file name, a GMTgrid or GMTimage object or a GDAL dataset

  • options: List of options. The accepted options are the ones of the gdal_translate utility. This list can be in the form of a vector of strings, or joined in a single string.

  • kwargs: Besides what was mentioned above one can also use meta=metadata, where metadata is a string vector with the form "NAME=...." for each of its elements. This data will be recognized by GDAL as Metadata. The kwargs may also contain the GMT region (-R), proj (-J), inc (-I) and save=fname options.

Returns

A GMT grid or Image, or a GDAL dataset (or nothing if file was writen on disk).

gdalwarp(indata, options=String[]; dest="/vsimem/tmp", kwargs...)

Image/Grid reprojection and warping function.

Parameters

  • indata: Input data. It can be a file name, a GMTgrid or GMTimage object or a GDAL dataset

  • options: List of options. The accepted options are the ones of the gdal_translate utility. This list can be in the form of a vector of strings, or joined in a single string. The accepted options are the ones of the gdalwarp utility.

  • kwargs: Besides what was mentioned above one can also use meta=metadata, where metadata is a string vector with the form "NAME=...." for each of its elements. This data will be recognized by GDAL as Metadata. The kwargs may also contain the GMT region (-R), proj (-J), inc (-I) and save=fname options.

Returns

A GMT grid or Image, or a GDAL dataset (or nothing if file was writen on disk).

gdaldem(dataset, method, options=String[]; dest="/vsimem/tmp", color=name|GMTcpt, kw...)

Tools to analyze and visualize DEMs.

Parameters

  • dataset The source dataset.

  • method the processing to apply (one of "hillshade", "slope", "aspect", "color-relief", "TRI", "TPI", "Roughness").

  • options List of options (potentially including filename and open options). The accepted options are the ones of the gdaldem utility.

Keyword Arguments

  • color color file (mandatory for "color-relief" processing, should be empty otherwise). If color is just a CPT name we compute a CPT from it and the input grid.

  • kw... keyword=value arguments when method is hillshade.

Returns

A GMT grid or Image, or a GDAL dataset (or nothing if file was writen on disk).

ogr2ogr(indata, options=String[]; dest="/vsimem/tmp", kwargs...)

Parameters

  • indata The source dataset.

  • options List of options (potentially including filename and open options). The accepted options are the ones of the gdalwarp utility.

  • kw are kwargs that may contain the GMT region (-R), proj (-J), inc (-I) and save=fname options

Returns

A GMT dataset, or a GDAL dataset (or nothing if file was writen on disk).

arccircle(x0, y0, radius, start_angle, end_angle; z0=0, inc=2, gdataset=false)

Parameters

  • x0: center X

  • y0: center Y

  • radius: radius of the circle.

  • start_angle: angle to first point on arc (clockwise of X-positive)

  • end_angle: angle to last point on arc (clockwise of X-positive)

Keywords

  • z0: center Z. Optional, if not provided the output is flat 2D

  • inc: the largest step in degrees along the arc. Default is 2 degree

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Returns

A GMT dataset or a GDAL IGeometry

arcellipse(x0, y0, primary_radius, secondary_radius, start_angle, end_angle; rotation=0, z0=0, inc=2, gdataset=false)

Parameters

  • x0: center X

  • y0: center Y

  • primary_radius: X radius of ellipse.

  • secondary_radius: Y radius of ellipse.

  • start_angle: angle to first point on arc (clockwise of X-positive)

  • end_angle: angle to last point on arc (clockwise of X-positive)

Keywords

  • rotation: rotation of the ellipse clockwise.

  • z0: center Z. Optional, if not provided the output is flat 2D

  • inc: the largest step in degrees along the arc. Default is 2 degree

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Returns

A GMT dataset or a GDAL IGeometry

boundary(geom; gdataset=false)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

A new geometry object is created and returned containing the boundary of the geometry on which the method is invoked.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

buffer(geom, dist::Real, quadsegs::Integer=30; gdataset=false)

Compute buffer of a geometry.

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • dist: the buffer distance to be applied. Should be expressed into the same unit as the coordinates of the geometry.

  • quadsegs: the number of segments used to approximate a 90 degree (quadrant) of curvature.

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix

Builds a new geometry containing the buffer region around the geometry on which it is invoked. The buffer is a polygon containing the region within the buffer distance of the original geometry.

Some buffer sections are properly described as curves, but are converted to approximate polygons. The quadsegs parameter can be used to control how many segments should be used to define a 90 degree curve - a quadrant of a circle. A value of 30 is a reasonable default. Large values result in large numbers of vertices in the resulting buffer geometry while small numbers reduce the accuracy of the result.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

buffergeo(D::GMTdataset; width=0, unit=:m, np=120, flatstart=false, flatend=false, epsg::Integer=0, tol=0.01)

or

buffergeo(line::Matrix; width=0, unit=:m, np=120, flatstart=false, flatend=false, proj::String="", epsg::Integer=0, tol=0.01)

or

buffergeo(fname::String; width=0, unit=:m, np=120, flatstart=false, flatend=false, proj::String="", epsg::Integer=0, tol=0.01)

Computes a buffer arround a poly-line. This calculation is performed on a ellipsoidal Earth (or other planet) using the GeographicLib (via PROJ4) so it should be very accurate.

Parameters

  • D | line | fname: - the geometry. This can either be a GMTdataset (or vector of it), a Mx2 matrix, the name of file that can be read as a GMTdataset by gmtread() or a GDAL AbstractDataset object

  • width: - the buffer width to be applied. Expressed meters (the default), km or Miles (see unit)

  • unit: - If width is not in meters use one of unit=:km, or unit=:Nautical or unit=:Miles

  • np: - Number of points into which circles are descretized (Default = 120)

  • flatstart - When computing buffers arround poly-lines make the start flat (no half-circle)

  • flatend - Same as flatstart but for the buffer end

  • proj - If line data is in Cartesians but with a known projection pass in a PROJ4 string to allow computing the buffer

  • epsg - Same as proj but using an EPSG code

  • tol - At the end simplify the buffer line with a Douglas-Peucker procedure. Use TOL=0 to NOT do the line simplification, or use any other value in degrees. Default computes it as 0.5% of buffer width.

Returns

A GMT dataset or a vector of it (when input is Vector{GMTdataset})

Example: Compute a buffer with 50000 m width

D = buffergeo([0 0; 10 10; 15 20], width=50000);
centroid(geom; gdataset=false)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of them), or a Matrix

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is a GMTdataset or Matrix

Compute the geometry centroid.

The centroid is not necessarily within the geometry.

(This method relates to the SFCOM ISurface::get_Centroid() method however the current implementation based on GEOS can operate on other geometry types such as multipoint, linestring, geometrycollection such as multipolygons. OGC SF SQL 1.1 defines the operation for surfaces (polygons). SQL/MM-Part 3 defines the operation for surfaces and multisurfaces (multipolygons).)

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

circgeo(lon, lat; radius=X, proj="", s_srs="", epsg=0, dataset=false, unit=:m, np=120, shape="")

or

circgeo(lonlat; radius=X, proj="", s_srs="", epsg=0, dataset=false, unit=:m, np=120, shape="")

Args:

  • lonlat: - longitude, latitude (degrees). If a Mx2 matrix, returns as many segments as number of rows. Use this to compute multiple shapes at different positions. In this case output type is

			always a vector of GMTdatasets.
  • radius: - The circle radius in meters (but see unit) or circumscribing circle for the other shapes

  • proj or s_srs: - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT

  • epsg: - Alternative way of specifying the projection [Default is WGS84]

  • dataset: - If true returns a GMTdataset instead of matrix (with single shapes)

  • unit: - If radius is not in meters use one of unit=:km, or unit=:Nautical or unit=:Miles

  • np: - Number of points into which the circle is descretized (Default = 120)

  • shape: - Optional string/symbol with "triangle", "square", "pentagon" or "hexagon" (or just the first char) to compute one of those geometries instead of a circle. np is ignored in these cases.

Returns

  • circ - A Mx2 matrix or GMTdataset with the circle coordinates

Example: Compute circle about the (0,0) point with a radius of 50 km

c = circgeo([0.,0], radius=50, unit=:k)
concavehull(geom, ratio, allow_holes::Bool=true; gdataset=false)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • ratio: Ratio of the area of the convex hull and the concave hull.

  • allow_holes: Whether holes are allowed.

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

contains(haystack::AbstractString, needle)

Return true if haystack contains needle. This is the same as occursin(needle, haystack), but is provided for consistency with startswith(haystack, needle) and endswith(haystack, needle).

See also occursin, in, issubset.

Examples

julia> contains("JuliaLang is pretty cool!", "Julia")
true

julia> contains("JuliaLang is pretty cool!", 'a')
true

julia> contains("aba", r"a.a")
true

julia> contains("abba", r"a.a")
false

Julia 1.5

The contains function requires at least Julia 1.5.

contains(needle)

Create a function that checks whether its argument contains needle, i.e. a function equivalent to haystack -> contains(haystack, needle).

The returned function is of type Base.Fix2{typeof(contains)}, which can be used to implement specialized methods.

contains(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if g1 contains g2.

convexhull(geom; gdataset=false)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

A new geometry object is created and returned containing the convex hull of the geometry on which the method is invoked.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

crosses(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if the geometries are crossing.

delaunay(geom::AbstractGeometry, tol::Real=0, onlyedges::Bool=true; gdataset=false)

Parameters

  • geom: the geometry.

  • tol: optional snapping tolerance to use for improved robustness

  • onlyedges: if true, will return a MULTILINESTRING, otherwise it will return a GEOMETRYCOLLECTION containing triangular POLYGONs.

  • gdataset: Returns a GDAL IGeometry even when input is GMTdataset or Matrix

Return a Delaunay triangulation of the vertices of the geometry.

difference(geom1, geom2; gdataset=false)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is GMTdataset or Matrix

Generates a new geometry which is the region of this geometry with the region of the other geometry removed.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

disjoint(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if the geometries are disjoint.

distance(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns the distance between the geometries or -1 if an error occurs.

dither(indata; n_colors=256, save="", gdataset=false)

Convert a 24bit RGB image to 8bit paletted.

  • Use the save=fname option to save the result to disk in a GeoTiff file fname. If fname has no extension a .tif one will be appended. Alternatively give file names with extension .png or .nc to save the file in one of those formats.

  • gdataset=true: to return a GDAL dataset. The default is to return a GMTimage object.

  • n_colors: Select the number of colors in the generated color table. Defaults to 256.

epsg2proj(code::Integer)

Convert a EPSG code into the PROJ4 form.

epsg2wkt(code::Integer, pretty::Bool=false)

Convert a EPSG code into the WKT form. Use pretty=true to return a more human readable text.

equals(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if the geometries are equivalent.

envelope(geom)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Computes and returns the bounding envelope for this geometry.

envelope3d(geom)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Computes and returns the bounding envelope (3D) for this geometry

GI = fillnodata(data::String; nodata=nothing, kwargs...) -> GMTgrid or GMTimage

Fill selected raster regions by interpolation from the edges.

Parameters

  • data: Input data. The file name of a grid or image that can be read with gmtread.

  • nodata: The nodata value that will be used to fill the regions. Otherwise use the nodata attribute of indata if it exists, or NaN if none of the above were set.

  • kwargs:

    • band: the band number. Default is first layer in indata

    • maxdist: the maximum number of cels to search in all directions to find values to interpolate from. Default, fills all.

    • nsmooth: the number of 3x3 smoothing filter passes to run (0 or more).

Returns

A GMTgrid or GMTimage object with the band nodata values filled by interpolation.

fillnodata!(data::GItype; nodata=nothing, kwargs...)

Fill selected raster regions by interpolation from the edges.

Parameters

  • data: Input data. It can be a file name, a GMTgrid or GMTimage object.

  • nodata: The nodata value that will be used to fill the regions. Otherwise use the nodata attribute of indata if it exists, or NaN if none of the above were set.

  • kwargs:

    • band: the band number. Default is first layer in indata

    • maxdist: the maximum number of cels to search in all directions to find values to interpolate from. Default, fills all.

    • nsmooth: the number of 3x3 smoothing filter passes to run (0 or more).

Returns

The modified input data

dest, azim = geod(lonlat, azim, distance; proj::String="", s_srs::String="", epsg::Integer=0, dataset=false, unit=:m)

Solve the direct geodesic problem.

Args:

  • lonlat: - longitude, latitude (degrees). This can be a vector or a matrix with one row only.

  • azimuth: - azimuth (degrees) ∈ [-540, 540)

  • distance: - distance to move from (lat,lon); can be vector and can be negative, Default is meters but see unit

  • proj or s_srs: - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT

  • epsg: - Alternative way of specifying the projection [Default is WGS84]

  • dataset: - If true returns a GMTdataset instead of matrix

  • unit: - If distance is not in meters use one of unit=:km, or unit=:Nautical or unit=:Miles

The distance argument can be a scalar, a Vector, a Vector{Vector} or an AbstractRange. The azimuth can be a scalar or a Vector.

When azimuth is a Vector we always return a GMTdataset with the multiple lines. Use this together with a non-scalar distance to get lines with multiple points along the line. The number of points along line does not need to be the same. For data, give the distance as a Vector{Vector} where each element of distance is a vector with the distances of the points along a line. In this case the number of distance elements must be equal to the number of azimuth.

Returns

  • dest - destination after moving for [distance] metres in [azimuth] direction.

  • azi - forward azimuth (degrees) at destination [dest].

Example: Compute two lines starting at (0,0) with lengths 111100 & 50000, heading at 15 and 45 degrees.

dest, = geod([0., 0], [15., 45], [[0, 10000, 50000, 111100.], [0., 50000]])[1]
function geodesic(D; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0, longest::Bool=false)

or

function geodesic(lon1, lat1, lon2, lat2; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0, longest::Bool=false)

Generate geodesic line(s) (shortest distace) on an ellipsoid. Input data can be two or more points. In later case each line segment is descretized at step increments,

Parameters

  • D: - the input points. This can either be a GMTdataset (or vector of it), a Mx2 matrix, the name of file that can be read as a GMTdataset by gmtread() or a GDAL AbstractDataset object

  • step: - Incremental distance at which the segment line is descretized in meters(the default), but see unit

  • unit: - If step is not in meters use one of unit=:km, or unit=:Nautical or unit=:Miles

  • np: - Number of intermediate points between poly-line vertices (alternative to step)

  • proj - If line data is in Cartesians but with a known projection pass in a PROJ4 string

  • epsg - Same as proj but using an EPSG code

  • longest - Compute the 'long way around' of the geodesic. That is, going from point A to B taking the longest path. But mind you that geodesics other than meridians and equator are not closed paths (see https://en.wikipedia.org/wiki/Geodesicsonanellipsoid). This line is obtained by computing the azimuth from A to B, at B, and start the line at B with that azimuth and go around the Earth till we reach, _close to A, but not exactly A (except for the simple meridian cases).

Returns

A Mx2 matrix with the lon lat of the points along the geodesic when input is a matrix. A GMT dataset or a vector of it (when input is Vector{GMTdataset}).

Example: Compute an geodesic between points (0,0) and (30,50) discretized at 100 km steps.

mat = geodesic([0 0; 30 50], step=100, unit=:k);
geodesicarea(geom)

Compute geometry area, considered as a surface on the underlying ellipsoid of the SRS attached to the geometry. The returned area will always be in square meters, and assumes that polygon edges describe geodesic lines on the ellipsoid. If the geometry' SRS is not a geographic one, geometries are reprojected to the underlying geographic SRS of the geometry' SRS.

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Returns the area of the geometry in square meters or a negative value in case of error for unsupported geometry types.

geomarea(geom)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Returns the area of the geometry or 0.0 for unsupported geometry types.

geomlength(geom)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

Returns the length of the geometry, or 0.0 for unsupported geometry types.

intersection(geom1, geom2; gdataset=false)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is GMTdataset or Matrix

Returns a new geometry representing the intersection of the geometries, or NULL if there is no intersection or an error occurs.

Generates a new geometry which is the region of intersection of the two geometries operated on. The intersects function can be used to test if two geometries intersect.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

intersects(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns whether the geometries intersect

Determines whether two geometries intersect. If GEOS is enabled, then this is done in rigorous fashion otherwise true is returned if the envelopes (bounding boxes) of the two geometries overlap.

dist, az1, az2 = invgeod(lonlat1::Vector{<:Real}, lonlat2::Vector{<:Real}; proj::String="",
                         s_srs::String="", epsg::Integer=0, backward=false)

Solve the inverse geodesic problem.

Args:

  • lonlat1: - coordinates of point 1 in the given projection (or a matrix with several points).

  • lonlat2: - coordinates of point 2 in the given projection (or a matrix with same size as lonlat1).

  • proj or s_srs: - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKT.

  • epsg: - Alternative way of specifying the projection [Default is WGS84].

  • backward: - If true, return backard azimuths.

Returns

dist - A scalar with the distance between point 1 and point 2 (meters). Or a vector when lonlat1|2 have more than one pair of points.

az1 - azimuth at point 1 (degrees) ∈ [-180, 180)

az2 - (forward) azimuth at point 2 (degrees) ∈ [-180, 180)

Remarks:

If either point is at a pole, the azimuth is defined by keeping the longitude fixed, writing lat = 90 +/- eps, and taking the limit as eps -> 0+.

lonlat2xy(lonlat::Matrix{<:Real}; t_srs, s_srs="+proj=longlat +datum=WGS84")

or

lonlat2xy(D::GMTdataset; t_srs, s_srs="+proj=longlat +datum=WGS84")

Computes the forward projection from LatLon to XY in the given projection. The input is assumed to be in WGS84. If it isn't, pass the appropriate projection info via the s_srs option (PROJ4, WKT, EPSG).

Parameters

  • lonlat: The input data. It can be a Matrix, or a GMTdataset (or vector of it)

  • t_srs: The destiny projection system. This can be a PROJ4, a WKT string or EPSG code

Returns

A Matrix if input is a Matrix or a GMTdadaset if input had that type

function loxodrome(lon1,lat1,lon2,lat2; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)

or

function loxodrome(D; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0)

Generate a loxodrome (rhumb line) on an ellipsoid. Input data can be two or more points. In later case each line segment is descretized at step increments,

Parameters

  • D: - the input points. This can either be a 2x2 matrix or a GMTdataset. Note that only the first 2 points are used.

  • step: - Incremental distance at which the segment line is descretized in meters(the default), but see unit

  • unit: - If step is not in meters use one of unit=:km, or unit=:Nautical or unit=:Miles

  • np: - Number of intermediate points to be generated between end points (alternative to step)

  • proj - If line data is in Cartesians but with a known projection pass in a PROJ4 string

  • epsg - Same as proj but using an EPSG code

Returns

A Mx2 matrix with the on lat of the points along the loxodrome when input is a matrix or the 2 pairs of points. A GMTdataset when the input is GMTdataset.

Example: Compute an loxodrome between points (0,0) and (30,50) discretized at 100 km steps.

loxo = loxodrome([0 0; 30 50], step=100, unit=:k);
loxodrome_direct(lon, lat, azimuth, distance, a=6378137.0, f=0.0033528106647474805)

Compute the direct problem of a loxodrome on the ellipsoid.

Given latitude and longitude of P1, azimuth a12 of the loxodrome P1P2 and the arc length s along the loxodrome curve, compute the latitude and longitude of P2.

Args:

  • lon, lat: - longitude, latitude (degrees) of starting point.

  • azimuth: - azimuth (degrees)

  • distance: - distance to move from (lat,lon) in meters

  • a - major axis of the ellipsoid (meters). Default values for WGS84

  • f - flattening od the ellipsoid (default = 1 / 298.257223563)

Returns

  • [lon lat] of destination after moving for [distance] metres in [azimuth] direction.

Example: Compute the end point at a bearing of 45 degrees 10000 meters from point 0,0

loxo = loxodrome_direct(0,0,45, 10000)
function loxodrome_inverse(lon1, lat1, lon2, lat2, a=6378137.0, f=0.0033528106647474805)

Compute the inverse problem of a loxodrome on the ellipsoid.

Given latitudes and longitudes of P1 and P2 on the ellipsoid, compute the azimuth a12 of the loxodrome P1P2, the arc length s along the loxodrome curve.

Args:

  • lon1, lat1, lon2, lat2: - longitude and latitude of starting and end points (degrees).

  • a - major axis of the ellipsoid (meters). Default values for WGS84

  • f - flattening od the ellipsoid (default = 1 / 298.257223563)

Returns

  • Distance (meters) and azimuth from P1 to P2

Example: Compute the distance and azimuth beyween points (0,0) and (5,5)

dist, azim = loxodrome_inverse(0,0,5,5)
O = gd2gmt(dataset; band=0, bands=[], sds=0, pad=0)

Convert a GDAL raster dataset into either a GMTgrid (if type is Int16 or Float) or a GMTimage type Use band to select a single band of the dataset. When you know that the dataset contains several bands of an image, use the kwarg bands with a vector the wished bands. By default it reads all bands of the image or grid object.

When DATASET is a string it may contain the file name or the name of a subdataset. In former case you can use the kwarg sds to selec the subdataset numerically. Alternatively, provide the full sds name. For files with sds with a scale_factor (e.g. MODIS data), that scale is applyied automaticaly.

Examples:

G = gd2gmt("AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc", sds=1);

or

G = gd2gmt("SUBDATASET1NAME=NETCDF:AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc:sst");

or

G = gd2gmt("NETCDF:AQUA_MODIS.20210228.L3m.DAY.NSST.sst.4km.NRT.nc:sst");

gdalshade(filename; kwargs...)

Create a shaded relief with the GDAL method (color image blended with shaded intensity).

  • kwargs hold the keyword=value to pass the arguments to gdaldem hillshade

Example: I = gdalshade("hawaii_south.grd", C="faa.cpt", zfactor=4);

Returns

A GMT RGB Image

gdalread(fname::AbstractString, opts=String[]; gdataset=false, kwargs...)

Read a raster or a vector file from a disk file and return the result either as a GMT type (the default) or a GDAL dataset.

  • fname: Input data. It can be a file name, a GMTgrid or GMTimage object or a GDAL dataset

  • opts: List of options. The accepted options are the ones of the gdal_translate utility. This list can be in the form of a vector of strings, or joined in a simgle string.

  • gdataset: If set to true forces the return of a GDAL dataset instead of a GMT type.

  • kwargs: This options accept the GMT region (-R) and increment (-I)

Returns

A GMT grid/image or a GDAL dataset

gdalwrite(fname::AbstractString, data, opts=String[]; kwargs...)

Write a raster or a vector file to disk

  • fname: Output file name. If not explicitly selected via opts the used driver will be picked from the file extension.

  • data: The data to be saved in file. It can be a GMT type or a GDAL dataset.

  • opts: List of options. The accepted options are the ones of the gdal_translate or ogr2ogr utility. This list can be in the form of a vector of strings, or joined in a simgle string.

  • kwargs: This options accept the GMT region (-R) and increment (-I)

or

gdalwrite(cube::GItype, fname::AbstractString, v=nothing; dim_name::String="time", dim_units::String="")

Write a MxNxP cube object to disk as a multilayered file.

  • cube: A GMTgrid or GMTimage cube

  • fname: The file name where to save the cube

  • v: A vector with the coordinates of the Z layers (if omitted create one as 1:size(cube,3))

  • dim_name: The name of the variable of the $vertical$ dimension.

  • dim_units: The units of the v vector. If not provided, use the cube.z_units if exist (GMTgrid only)

ds = gmt2gd(GI)

Create GDAL dataset from the contents of GI that can be either a Grid or an Image

ds = gmt2gd(D, save="", geometry="")

Create GDAL dataset from the contents of D, which can be a GMTdataset, a vector of GMTdataset ir a MxN array. The save keyword instructs GDAL to save the contents as an OGR file. Format is determined by file estension. geometry can be a string with "polygon", where file will be converted to polygon/multipolygon depending on D is a single or a multi-segment object, or "point" to convert to a multipoint geometry.

function geodesic(D; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0, longest::Bool=false)

or

function geodesic(lon1, lat1, lon2, lat2; step=0, unit=:m, np=0, proj::String="", epsg::Integer=0, longest::Bool=false)

Generate geodesic line(s) (shortest distace) on an ellipsoid. Input data can be two or more points. In later case each line segment is descretized at step increments,

Parameters

  • D: - the input points. This can either be a GMTdataset (or vector of it), a Mx2 matrix, the name of file that can be read as a GMTdataset by gmtread() or a GDAL AbstractDataset object

  • step: - Incremental distance at which the segment line is descretized in meters(the default), but see unit

  • unit: - If step is not in meters use one of unit=:km, or unit=:Nautical or unit=:Miles

  • np: - Number of intermediate points between poly-line vertices (alternative to step)

  • proj - If line data is in Cartesians but with a known projection pass in a PROJ4 string

  • epsg - Same as proj but using an EPSG code

  • longest - Compute the 'long way around' of the geodesic. That is, going from point A to B taking the longest path. But mind you that geodesics other than meridians and equator are not closed paths (see https://en.wikipedia.org/wiki/Geodesicsonanellipsoid). This line is obtained by computing the azimuth from A to B, at B, and start the line at B with that azimuth and go around the Earth till we reach, _close to A, but not exactly A (except for the simple meridian cases).

Returns

A Mx2 matrix with the lon lat of the points along the geodesic when input is a matrix. A GMT dataset or a vector of it (when input is Vector{GMTdataset}).

Example: Compute an geodesic between points (0,0) and (30,50) discretized at 100 km steps.

mat = geodesic([0 0; 30 50], step=100, unit=:k);
overlaps(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if the geometries overlap.

pointalongline(geom, distance::Real; gdataset=false)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • distance: distance along the curve at which to sample position. This distance should be between zero and geomlength() for this curve.

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Fetch point (or NULL) at given distance along curve.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

polygonize(geom; gdataset=false)

Parameters

  • geom: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • gdataset: Returns a GDAL IGeometry even when input are GMTdataset or Matrix

Polygonizes a set of sparse edges.

A new geometry object is created and returned containing a collection of reassembled Polygons: NULL will be returned if the input collection doesn't correspond to a MultiLinestring, or when reassembling Edges into Polygons is impossible due to topological inconsistencies.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

D = polygonize(data::GItype; kwargs...) -> Vector{GMTdataset}

This method, which uses the GDAL GDALPolygonize function, creates vector polygons for all connected regions of pixels in the raster sharing a common pixel/cell value. The input may be a grid or an image. This function can be rather slow as it picks lots of polygons in pixels with slightly different values at transitions between colors. Its natural use is to digitize masks images.

Parameters

  • data: Input data. It can be a GMTgrid or GMTimage object.

  • kwargs:

    • min, nmin, npixels or ncells: The minimum number of cells/pixels for a polygon to be retained. Default is 1. This can be set to filter out small polygons.

    • min_area: Minimum area in m2 for a polygon to be retained. This option takes precedence over the one above that is based in the counting of cells. Note also that this is an approximate value because at this point we still don't know exactly the latitudes of the polygons.

    • max_area: Maximum area in m2 for a polygon to be retained.

    • simplify: Apply the Douglas-Peucker line simplification algorithm to the poligons. Provide a tolerence in meters. For example: simplify=0.5. But be warned that this is a risky option since a too large tolerance

can lead to loss of otherwise good polygons. A good rule of thumb is to use the cell size for the tolerance.
And in fact that is what we do when using `simplify=:auto`.
  • sort: If true, will sort polygons by pixel count. Default is the order that GDAL decides internally.

polyunion(geom1, geom2; gdataset=false)

Computes a new geometry representing the union of the geometries.

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is GMTdataset or Matrix

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

proj2wkt(proj4_str::String, pretty::Bool=false)

Convert a PROJ4 string into the WKT form. Use pretty=true to return a more human readable text.

setproj!(type, proj)

Set a referencing system to the type object. This object can be a GMTgrid, a GMTimage, a GMTdataset or an AbstractDataset.

  • proj Is either a Proj4 string or a WKT. Alternatively, it can also be another grid, image or dataset type, in which case its referencing system is copied into type

simplify(geom::AbstractGeometry, tol::Real; gdataset=false)

Compute a simplified geometry.

Parameters

  • geom: the geometry.

  • tol: the distance tolerance for the simplification.

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is GMTdataset or Matrix

symdifference(geom1, geom2; gdataset=false)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Keywords

  • gdataset: Returns a GDAL IGeometry even when input is GMTdataset or Matrix

Generates a new geometry representing the symmetric difference of the geometries or NULL if the difference is empty or an error occurs.

Returns

A GMT dataset when input is a Matrix or a GMT type (except if gdaset=true), or a GDAL IGeometry otherwise

touches(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if the geometries are touching.

within(geom1, geom2)

Parameters

  • geom1: the geometry. This can either be a GDAL AbstractGeometry or a GMTdataset (or vector of it), or a Matrix

  • geom2: Second geometry. AbstractGeometry if geom1::AbstractGeometry or a GMTdataset/Matrix if geom1 is GMT type

Returns true if g1 is contained within g2.

wkt2proj(wkt_str::String)

Convert a WKT SRS string into the PROJ4 form.

xy2lonlat(xy::Matrix{<:Real}, s_srs=""; s_srs="", t_srs="+proj=longlat +datum=WGS84")

or

xy2lonlat(D::GMTdataset, s_srs=""; s_srs="", t_srs="+proj=longlat +datum=WGS84")

Computes the inverse projection from XY to LonLat in the given projection. The output is assumed to be in WGS84. If that isn't right, pass the appropriate projection info via the t_srs option (PROJ4, WKT, EPSG).

Parameters

  • xy: The input data. It can be a Matrix, or a GMTdataset (or vector of it)

  • s_srs: The data projection system. This can be a PROJ4, a WKT string or EPSG code

  • t_srs: The target SRS. If the default is not satisfactory, provide a new projection info (PROJ4, WKT, EPSG)

Returns

A Matrix if input is a Matrix or a GMTdadaset if input had that type