# 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 simgle 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.

### 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 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 (potentially including filename and open options). The accepted options are the ones of the gdalwarp utility.`kwargs`

: Are kwargs that may 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", colorfile=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

`colorfile`

color file (mandatory for "color-relief" processing, should be empty otherwise).`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)`

### 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.`gdataset`

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

Compute buffer of geometry.

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 it), or a Matrix`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

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

`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

`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`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

`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

`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/Geodesics*on*an*ellipsoid). 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);`

`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`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

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

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]

### 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("SUBDATASET_1_NAME=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(data, fname::AbstractString, opts=String[]; kwargs...)`

or

`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, 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/Geodesics*on*an*ellipsoid). 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

`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

`gdaset=true`

), or a GDAL IGeometry otherwise

`polyunion(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`gdataset`

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

Computes a new geometry representing the union of the geometries.

### Returns

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

`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

These docs were autogenerated using GMT: v0.44.4