Skip to content

Commit e204dda

Browse files
committed
Merge branch 'main' into wrap/gmtselect
2 parents 4d003de + 70f6f09 commit e204dda

File tree

11 files changed

+435
-4
lines changed

11 files changed

+435
-4
lines changed

.github/workflows/publish-to-pypi.yml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@ jobs:
3131
- name: Set up Python
3232
uses: actions/setup-python@v2.2.2
3333
with:
34-
python-version: 3.9
34+
python-version: '3.10'
3535

3636
- name: Install dependencies
3737
run: python -m pip install setuptools wheel

.github/workflows/style_checks.yaml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ jobs:
2222
- name: Set up Python
2323
uses: actions/setup-python@v2.2.2
2424
with:
25-
python-version: 3.9
25+
python-version: '3.10'
2626

2727
- name: Install packages
2828
run: |

doc/api/index.rst

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -83,6 +83,7 @@ Operations on tabular data:
8383
blockmedian
8484
blockmode
8585
nearneighbor
86+
project
8687
select
8788
sph2grd
8889
sphdistance

pygmt/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@
5050
info,
5151
makecpt,
5252
nearneighbor,
53+
project,
5354
select,
5455
sph2grd,
5556
sphdistance,

pygmt/accessors.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@ def __init__(self, xarray_obj):
3434
self._registration, self._gtype = map(
3535
int, grdinfo(self._source, per_column="n", o="10,11").split()
3636
)
37-
except KeyError:
37+
except (KeyError, ValueError):
3838
self._registration = 0 # Default to Gridline registration
3939
self._gtype = 0 # Default to Cartesian grid type
4040

pygmt/src/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
from pygmt.src.nearneighbor import nearneighbor
3737
from pygmt.src.plot import plot
3838
from pygmt.src.plot3d import plot3d
39+
from pygmt.src.project import project
3940
from pygmt.src.rose import rose
4041
from pygmt.src.select import select
4142
from pygmt.src.solar import solar

pygmt/src/grd2xyz.py

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,9 +18,17 @@
1818

1919
@fmt_docstring
2020
@use_alias(
21+
C="cstyle",
2122
R="region",
2223
V="verbose",
24+
W="weight",
25+
Z="convention",
26+
b="binary",
27+
d="nodata",
28+
f="coltypes",
29+
h="header",
2330
o="outcols",
31+
s="skiprows",
2432
)
2533
@kwargs_to_strings(R="sequence", o="sequence_comma")
2634
def grd2xyz(grid, output_type="pandas", outfile=None, **kwargs):
@@ -48,12 +56,65 @@ def grd2xyz(grid, output_type="pandas", outfile=None, **kwargs):
4856
- ``file`` - ASCII file (requires ``outfile``)
4957
outfile : str
5058
The file name for the output ASCII file.
59+
cstyle : str
60+
[**f**\|\ **i**].
61+
Replace the x- and y-coordinates on output with the corresponding
62+
column and row numbers. These start at 0 (C-style counting); append
63+
**f** to start at 1 (Fortran-style counting). Alternatively, append
64+
**i** to write just the two columns *index* and *z*, where *index*
65+
is the 1-D indexing that GMT uses when referring to grid nodes.
5166
{R}
5267
Adding ``region`` will select a subsection of the grid. If this
5368
subsection exceeds the boundaries of the grid, only the common region
5469
will be output.
70+
weight : str
71+
[**a**\ [**+u**\ *unit*]\|\ *weight*].
72+
Write out *x,y,z,w*\ , where *w* is the supplied *weight* (or 1 if not
73+
supplied) [Default writes *x,y,z* only]. Choose **a** to compute
74+
weights equal to the area each node represents. For Cartesian grids
75+
this is simply the product of the *x* and *y* increments (except for
76+
gridline-registered grids at all sides [half] and corners [quarter]).
77+
For geographic grids we default to a length unit of **k**. Change
78+
this by appending **+u**\ *unit*. For such grids, the area
79+
varies with latitude and also sees special cases for
80+
gridline-registered layouts at sides, corners, and poles.
5581
{V}
82+
convention : str
83+
[*flags*].
84+
Write a 1-column ASCII [or binary] table. Output will be organized
85+
according to the specified ordering convention contained in *flags*.
86+
If data should be written by rows, make *flags* start with
87+
**T** (op) if first row is y = ymax or
88+
**B** (ottom) if first row is y = ymin. Then,
89+
append **L** or **R** to indicate that first element should start at
90+
left or right end of row. Likewise for column formats: start with
91+
**L** or **R** to position first column, and then append **T** or
92+
**B** to position first element in a row. For gridline registered
93+
grids: If grid is periodic in x but the written data should not
94+
contain the (redundant) column at x = xmax, append **x**. For grid
95+
periodic in y, skip writing the redundant row at y = ymax by
96+
appending **y**. If the byte-order needs to be swapped, append
97+
**w**. Select one of several data types (all binary except **a**):
98+
99+
* **a** ASCII representation of a single item per record
100+
* **c** int8_t, signed 1-byte character
101+
* **u** uint8_t, unsigned 1-byte character
102+
* **h** int16_t, short 2-byte integer
103+
* **H** uint16_t, unsigned short 2-byte integer
104+
* **i** int32_t, 4-byte integer
105+
* **I** uint32_t, unsigned 4-byte integer
106+
* **l** int64_t, long (8-byte) integer
107+
* **L** uint64_t, unsigned long (8-byte) integer
108+
* **f** 4-byte floating point single precision
109+
* **d** 8-byte floating point double precision
110+
111+
Default format is scanline orientation of ASCII numbers: **TLa**.
112+
{b}
113+
{d}
114+
{f}
115+
{h}
56116
{o}
117+
{s}
57118
58119
Returns
59120
-------

pygmt/src/project.py

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
"""
2+
project - Project data onto lines or great circles, or generate tracks.
3+
"""
4+
import pandas as pd
5+
from pygmt.clib import Session
6+
from pygmt.exceptions import GMTInvalidInput
7+
from pygmt.helpers import (
8+
GMTTempFile,
9+
build_arg_string,
10+
fmt_docstring,
11+
kwargs_to_strings,
12+
use_alias,
13+
)
14+
15+
16+
@fmt_docstring
17+
@use_alias(
18+
A="azimuth",
19+
C="center",
20+
E="endpoint",
21+
F="convention",
22+
G="generate",
23+
L="length",
24+
N="flat_earth",
25+
Q="unit",
26+
S="sort",
27+
T="pole",
28+
V="verbose",
29+
W="width",
30+
Z="ellipse",
31+
f="coltypes",
32+
)
33+
@kwargs_to_strings(E="sequence", L="sequence", T="sequence", W="sequence", C="sequence")
34+
def project(data=None, x=None, y=None, z=None, outfile=None, **kwargs):
35+
r"""
36+
Project data onto lines or great circles, or generate tracks.
37+
38+
Project reads arbitrary :math:`(x, y [, z])` data and returns any
39+
combination of :math:`(x, y, z, p, q, r, s)`, where :math:`(p, q)` are the
40+
coordinates in the projection, :math:`(r, s)` is the position in the
41+
:math:`(x, y)` coordinate system of the point on the profile (:math:`q = 0`
42+
path) closest to :math:`(x, y)`, and :math:`z` is all remaining columns in
43+
the input (beyond the required :math:`x` and :math:`y` columns).
44+
45+
Alternatively, ``project`` may be used to generate
46+
:math:`(r, s, p)` triples at equal increments along a profile using the
47+
``generate`` parameter. In this case, the value of ``data`` is ignored
48+
(you can use, e.g., ``data=None``).
49+
50+
Projections are defined in any (but only) one of three ways:
51+
52+
1. By a ``center`` and an ``azimuth`` in degrees clockwise from North.
53+
2. By a ``center`` and ``endpoint`` of the projection path.
54+
3. By a ``center`` and a ``pole`` position.
55+
56+
To spherically project data along a great circle path, an oblique
57+
coordinate system is created which has its equator along that path, and the
58+
zero meridian through the Center. Then the oblique longitude (:math:`p`)
59+
corresponds to the distance from the Center along the great circle, and the
60+
oblique latitude (:math:`q`) corresponds to the distance perpendicular to
61+
the great circle path. When moving in the increasing (:math:`p`) direction,
62+
(toward B or in the azimuth direction), the positive (:math:`q`) direction
63+
is to your left. If a Pole has been specified, then the positive
64+
(:math:`q`) direction is toward the pole.
65+
66+
To specify an oblique projection, use the ``pole`` option to set
67+
the pole. Then the equator of the projection is already determined and the
68+
``center`` option is used to locate the :math:`p = 0` meridian. The center
69+
*cx/cy* will be taken as a point through which the :math:`p = 0` meridian
70+
passes. If you do not care to choose a particular point, use the South pole
71+
(*cx* = 0, *cy* = -90).
72+
73+
Data can be selectively windowed by using the ``length`` and ``width``
74+
options. If ``width`` is used, the projection width is set to use only
75+
data with :math:`w_{{min}} < q < w_{{max}}`. If ``length`` is set, then
76+
the length is set to use only those data with
77+
:math:`l_{{min}} < p < l_{{max}}`. If the ``endpoint`` option
78+
has been used to define the projection, then ``length="w"`` may be used to
79+
window the length of the projection to exactly the span from O to B.
80+
81+
Flat Earth (Cartesian) coordinate transformations can also be made. Set
82+
``flat_earth=True`` and remember that azimuth is clockwise from North (the
83+
y axis), NOT the usual cartesian theta, which is counterclockwise from the
84+
x axis. azimuth = 90 - theta.
85+
86+
No assumptions are made regarding the units for
87+
:math:`x, y, r, s, p, q, dist, l_{{min}}, l_{{max}}, w_{{min}}, w_{{max}}`.
88+
If -Q is selected, map units are assumed and :math:`x, y, r, s` must be in
89+
degrees and :math:`p, q, dist, l_{{min}}, l_{{max}}, w_{{min}}, w_{{max}}`
90+
will be in km.
91+
92+
Calculations of specific great-circle and geodesic distances or for
93+
back-azimuths or azimuths are better done using :gmt-docs:`mapproject` as
94+
project is strictly spherical.
95+
96+
{aliases}
97+
98+
Parameters
99+
----------
100+
data : str or {table-like}
101+
Pass in (x, y, z) or (longitude, latitude, elevation) values by
102+
providing a file name to an ASCII data table, a 2D
103+
{table-classes}.
104+
105+
center : str or list
106+
*cx*/*cy*.
107+
Set the origin of the projection, in Definition 1 or 2. If
108+
Definition 3 is used, then *cx/cy* are the coordinates of a
109+
point through which the oblique zero meridian (:math:`p = 0`) should
110+
pass. The *cx/cy* is not required to be 90 degrees from the pole.
111+
112+
azimuth : float or str
113+
Define the azimuth of the projection (Definition 1).
114+
115+
endpoint : str or list
116+
*bx*/*by*.
117+
Define the end point of the projection path (Definition 2).
118+
119+
convention : str
120+
Specify the desired output using any combination of **xyzpqrs**, in
121+
any order [Default is **xypqrsz**]. Do not space between the letters.
122+
Use lower case. The output will be columns of values corresponding to
123+
your ``convention``. The **z** flag is special and refers to all
124+
numerical columns beyond the leading **x** and **y** in your input
125+
record. The **z** flag also includes any trailing text (which is
126+
placed at the end of the record regardless of the order of **z** in
127+
``convention``). **Note**: If ``generate`` is True, then the output
128+
order is hardwired to be **rsp** and ``convention`` is not allowed.
129+
130+
generate : str
131+
*dist* [/*colat*][**+c**\|\ **h**].
132+
Create :math:`(r, s, p)` output data every *dist* units of :math:`p`
133+
(See `unit` option). Alternatively, append */colat* for a small
134+
circle instead [Default is a colatitude of 90, i.e., a great circle].
135+
If setting a pole with ``pole`` and you want the small circle to go
136+
through *cx*/*cy*, append **+c** to compute the required colatitude.
137+
Use ``center`` and ``endpoint`` to generate a circle that goes
138+
through the center and end point. Note, in this case the center and
139+
end point cannot be farther apart than :math:`2|\mbox{{colat}}|`.
140+
Finally, if you append **+h** then we will report the position of
141+
the pole as part of the segment header [Default is no header].
142+
Note: No input is read and the value of ``data``, ``x``, ``y``,
143+
and ``z`` is ignored if ``generate`` is used.
144+
145+
length : str or list
146+
[**w**\|\ *l_min*/*l_max*].
147+
Project only those data whose *p* coordinate is
148+
within :math:`l_{{min}} < p < l_{{max}}`. If ``endpoint`` has been set,
149+
then you may alternatively use **w** to stay within the distance from
150+
``center`` to ``endpoint``.
151+
152+
flat_earth : bool
153+
Make a Cartesian coordinate transformation in the plane.
154+
[Default is ``False``; plane created with spherical trigonometry.]
155+
156+
unit : bool
157+
Set units for :math:`x, y, r, s` degrees and
158+
:math:`p, q, dist, l_{{min}}, l_{{max}}, w_{{min}}, {{w_max}}` to km.
159+
[Default is ``False``; all arguments use the same units]
160+
161+
sort : bool
162+
Sort the output into increasing :math:`p` order. Useful when projecting
163+
random data into a sequential profile.
164+
165+
pole : str or list
166+
*px*/*py*.
167+
Set the position of the rotation pole of the projection.
168+
(Definition 3).
169+
170+
{V}
171+
172+
width : str or list
173+
*w_min*/*w_max*.
174+
Project only those data whose :math:`q` coordinate is
175+
within :math:`w_{{min}} < q < w_{{max}}`.
176+
177+
ellipse : str
178+
*major*/*minor*/*azimuth* [**+e**\|\ **n**].
179+
Used in conjunction with ``center`` (sets its center) and ``generate``
180+
(sets the distance increment) to create the coordinates of an ellipse
181+
with *major* and *minor* axes given in km (unless ``flat_earth`` is
182+
given for a Cartesian ellipse) and the *azimuth* of the major axis in
183+
degrees. Append **+e** to adjust the increment set via ``generate`` so
184+
that the the ellipse has equal distance increments [Default uses the
185+
given increment and closes the ellipse]. Instead, append **+n** to set
186+
a specific number of unique equidistant data via ``generate``. For
187+
degenerate ellipses you can just supply a single *diameter* instead. A
188+
geographic diameter may be specified in any desired unit other than km
189+
by appending the unit (e.g., 3d for degrees) [Default is km];
190+
the increment is assumed to be in the same unit. **Note**:
191+
For the Cartesian ellipse (which requires ``flat_earth``), the
192+
*direction* is counter-clockwise from the horizontal instead of an
193+
*azimuth*.
194+
195+
outfile : str
196+
The file name for the output ASCII file.
197+
198+
{f}
199+
200+
Returns
201+
-------
202+
track: pandas.DataFrame or None
203+
Return type depends on whether the ``outfile`` parameter is set:
204+
205+
- :class:`pandas.DataFrame` table with (x, y, ..., newcolname) if
206+
``outfile`` is not set
207+
- None if ``outfile`` is set (output will be stored in file set
208+
by ``outfile``)
209+
"""
210+
211+
if "C" not in kwargs:
212+
raise GMTInvalidInput("The `center` parameter must be specified.")
213+
if "G" not in kwargs and data is None:
214+
raise GMTInvalidInput(
215+
"The `data` parameter must be specified unless `generate` is used."
216+
)
217+
if "G" in kwargs and "F" in kwargs:
218+
raise GMTInvalidInput(
219+
"The `convention` parameter is not allowed with `generate`."
220+
)
221+
222+
with GMTTempFile(suffix=".csv") as tmpfile:
223+
if outfile is None: # Output to tmpfile if outfile is not set
224+
outfile = tmpfile.name
225+
with Session() as lib:
226+
if "G" not in kwargs:
227+
# Choose how data will be passed into the module
228+
table_context = lib.virtualfile_from_data(
229+
check_kind="vector", data=data, x=x, y=y, z=z, required_z=False
230+
)
231+
232+
# Run project on the temporary (csv) data table
233+
with table_context as infile:
234+
arg_str = " ".join(
235+
[infile, build_arg_string(kwargs), "->" + outfile]
236+
)
237+
else:
238+
arg_str = " ".join([build_arg_string(kwargs), "->" + outfile])
239+
lib.call_module(module="project", args=arg_str)
240+
241+
# if user did not set outfile, return pd.DataFrame
242+
if outfile == tmpfile.name:
243+
if "G" in kwargs:
244+
column_names = list("rsp")
245+
result = pd.read_csv(tmpfile.name, sep="\t", names=column_names)
246+
else:
247+
result = pd.read_csv(tmpfile.name, sep="\t", header=None, comment=">")
248+
# return None if outfile set, output in outfile
249+
elif outfile != tmpfile.name:
250+
result = None
251+
252+
return result

0 commit comments

Comments
 (0)