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 datasetoptions
: 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 usemeta=metadata
, wheremetadata
is a string vector with the form "NAME=...." for each of its elements. This data will be recognized by GDAL as Metadata. Thekwargs
may also contain the GMT region (-R), proj (-J), inc (-I) andsave=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 datasetoptions
: 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 usemeta=metadata
, wheremetadata
is a string vector with the form "NAME=...." for each of its elements. This data will be recognized by GDAL as Metadata. Thekwargs
may also contain the GMT region (-R), proj (-J), inc (-I) andsave=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). Ifcolor
is just a CPT name we compute a CPT from it and the input grid.kw...
keyword=value arguments whenmethod
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) andsave=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 Xy0
: center Yradius
: 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 2Dinc
: the largest step in degrees along the arc. Default is 2 degreegdataset
: 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 Xy0
: center Yprimary_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 2Dinc
: the largest step in degrees along the arc. Default is 2 degreegdataset
: 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 Matrixgdataset
: 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 Matrixdist
: 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 bygmtread()
or a GDAL AbstractDataset objectwidth
: - the buffer width to be applied. Expressed meters (the default), km or Miles (seeunit
)unit
: - Ifwidth
is not in meters use one ofunit=:km
, orunit=:Nautical
orunit=: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 asflatstart
but for the buffer endproj
- If line data is in Cartesians but with a known projection pass in a PROJ4 string to allow computing the bufferepsg
- Same asproj
but using an EPSG codetol
- 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 seeunit
) or circumscribing circle for the other shapesproj
ors_srs
: - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKTepsg
: - Alternative way of specifying the projection [Default is WGS84]dataset
: - If true returns a GMTdataset instead of matrix (with single shapes)unit
: - Ifradius
is not in meters use one ofunit=:km
, orunit=:Nautical
orunit=: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 Matrixratio
: 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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixgdataset
: 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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 robustnessonlyedges
: iftrue
, 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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 filefname
. Iffname
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 withgmtread
.nodata
: The nodata value that will be used to fill the regions. Otherwise use thenodata
attribute ofindata
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 thenodata
attribute ofindata
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 seeunit
proj
ors_srs
: - the given projection whose ellipsoid we move along. Can be a proj4 string or an WKTepsg
: - Alternative way of specifying the projection [Default is WGS84]dataset
: - If true returns a GMTdataset instead of matrixunit
: - Ifdistance
is not in meters use one ofunit=:km
, orunit=:Nautical
orunit=: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 bygmtread()
or a GDAL AbstractDataset objectstep
: - Incremental distance at which the segment line is descretized in meters(the default), but seeunit
unit
: - Ifstep
is not in meters use one ofunit=:km
, orunit=:Nautical
orunit=:Miles
np
: - Number of intermediate points between poly-line vertices (alternative tostep
)proj
- If line data is in Cartesians but with a known projection pass in a PROJ4 stringepsg
- Same asproj
but using an EPSG codelongest
- 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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 aslonlat1
).proj
ors_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
: - Iftrue
, 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 seeunit
unit
: - Ifstep
is not in meters use one ofunit=:km
, orunit=:Nautical
orunit=:Miles
np
: - Number of intermediate points to be generated between end points (alternative tostep
)proj
- If line data is in Cartesians but with a known projection pass in a PROJ4 stringepsg
- Same asproj
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 metersa
- major axis of the ellipsoid (meters). Default values for WGS84f
- 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 WGS84f
- 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 datasetopts
: 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 totrue
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 viaopts
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 cubefname
: The file name where to save thecube
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 thev
vector. If not provided, use thecube.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 bygmtread()
or a GDAL AbstractDataset objectstep
: - Incremental distance at which the segment line is descretized in meters(the default), but seeunit
unit
: - Ifstep
is not in meters use one ofunit=:km
, orunit=:Nautical
orunit=:Miles
np
: - Number of intermediate points between poly-line vertices (alternative tostep
)proj
- If line data is in Cartesians but with a known projection pass in a PROJ4 stringepsg
- Same asproj
but using an EPSG codelongest
- 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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixdistance
: 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 Matrixgdataset
: 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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 intotype
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 Matrixgeom2
: Second geometry. AbstractGeometry ifgeom1::AbstractGeometry
or a GMTdataset/Matrix ifgeom1
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 codet_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: v1.11.0