Support various datetime types as input

[seisman]

Description of proposed changes

This PR allows PyGMT accept vectors of various datetime types so that PyGMT can plot datetime axis.

Vectors of following datetime types are supported:

• numpy.datetime64
• pandas.DateTimeIndex
• Raw datetime strings in both ISO and non-ISO formats, e.g., 2010-01-01, 1/1/2018, Jul 5, 2019
• Python's built-in datetime and date

TODO:

• Wait for the GMT 6.1.0 release
• Add more tests
• Add a tutorial for calendar axis
• Test for xarray.DataArray

Known issues but won't be fixed in this PR:

• The region argument doesn't work with numpy datetime64 objects, have to convert to a string using np.datetime_as_string in order to set the map frame bounds.
• Support timedelta64, e.g. GPS time from a certain epoch.

Fixes #242.

Example script:

import datetime
import pygmt
import numpy as np
import pandas as pd

fig = pygmt.Figure()
fig.basemap(projection="X15c/5c", region="2010-01-01/2020-01-01/0/10", frame=True)

# numpy.datetime64 types
x = np.array(
    ["2010-06-01", "2011-06-01T12", "2012-01-01T12:34:56"], dtype="datetime64"
)
y = [1.0, 2.0, 3.0]
fig.plot(x, y, style="c0.2c", pen="1p")

# pandas.DatetimeIndex
x = pd.date_range("2013", freq="YS", periods=3)
y = [4, 5, 6]
fig.plot(x, y, style="t0.2c", pen="1p")

# xarray.DataArray
x = xr.DataArray(data=pd.date_range(start="2015-03", freq="QS", periods=3))
y = [7.5, 6, 4.5]
fig.plot(x, y, style="s0.2c", pen="1p")

# raw datetime strings
x = ["2016-02-01", "2017-03-04T00:00"]
y = [7, 8]
fig.plot(x, y, style="a0.2c", pen="1p")

# the Python built-in datetime and date
x = [datetime.date(2018, 1, 1), datetime.datetime(2019, 1, 1)]
y = [8.5, 9.5]
fig.plot(x, y, style="i0.2c", pen="1p")

fig.savefig("datetime.pdf")

Output:

Reminders

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst.
• Write detailed docstrings for all functions/methods.
• If adding new functionality, add an example to docstrings or tutorials.

[seisman]

@GenericMappingTools/python @GenericMappingTools/python-contributors This PR is almost done. If you have the GMT latest master branch installed, please try this branch and leave your comments.

[weiji14]

Some instructions on how to test this branch:

curl https://raw.githubusercontent.com/GenericMappingTools/gmt/master/ci/build-gmt-master.sh | bash
conda activate --name pygmt-test  # or some virtual env you have
pip install https://github.com/GenericMappingTools/pygmt/archive/datetime-input.zip

and then in python:

import os

os.environ["GMT_LIBRARY_PATH"] = "/home/username/gmt-install-dir/lib"

import pygmt

pygmt.print_clib_info()  # check that you're using GMT from master

and then try to test out plot with datetimes! I'll have a go at this in a bit with some 'real' data.

[weiji14]

Ok, I tried installing it in a different conda environment and it works now! Below is my convoluted real-data example from Antarctica (longitude=46.522298201524514, latitude=-73.05774246536473).

Code:

# %%
import os

os.environ["GMT_LIBRARY_PATH"] = "/home/username/gmt-install-dir/lib"

import numpy as np
import pygmt

pygmt.show_versions()

# %%
x = np.array(
    [
        np.datetime64(dt)
        for dt in [
            "NaT",
            "2019-01-28T07:51:16.582496785",
            "2019-04-29T03:30:57.620034382",
            "2019-07-28T23:10:36.096639410",
            "2019-10-27T18:50:32.615134843",
            "2020-01-26T14:30:18.835180975",
        ]
    ]
)
y = np.array(
    [np.NaN, 3126.60298909, 3126.67885045, 3126.67276984, 3126.7217197, 3126.68279204,]
)

# %%
# Mask out NaT values
mask = ~np.isnan(y)
x = x[mask]
y = y[mask]

# %%
dates = np.datetime_as_string(x, unit="ns")
region = [dates[0], dates[-1], np.nanmin(y) - 0.1, np.nanmax(y) + 0.1]

# %%
fig = pygmt.Figure()
fig.basemap(
    projection="X15c/5c",
    region=region,
    frame=["WSne", "xaf+lDateTime", "yaf+lHeight(m)"],
)
fig.plot(x=x, y=y, style="t1c", pen="1p")
fig.savefig("height_over_time.png")
fig.show()

Glad to see that these work off the shelf:

• Passing in xarray.DataArrays to plot works (thanks to NEP18's __array__ functionality).
• ISO datetimes down to nanoseconds are supported by GMT!

Places to improve:

• The region argument doesn't work with numpy datetime64 objects, have to convert to a string using np.datetime_as_string in order to set the map frame bounds.
• Having even just one NaT value will raise an error, as mentioned in Support datetime data types as input #242 (comment), should be a warning ideally, but it's an upstream issue. Can be tested by commenting out the 'mask' code above.
• Support timedelta64, e.g. GPS time from a certain epoch. More of a nice-to-have really, could be done in a separate PR once this datetime one is merged in.

[seisman]

• Passing in xarray.DataArrays to plot works (thanks to NEP18's __array__ functionality).

I'll add it to the test

• ISO datetimes down to nanoseconds are supported by GMT!

Perhaps not true. GMT converts ISO datetime to double values internally. 2019-01-28 07:51:16.582496785 is stored as a double value 1548661876.582496881484985 (printed with %30.15lf). So GMT only supports up to microsecond, as limited by the precision of double type. Maybe I'm wrong here.

• The region argument doesn't work with numpy datetime64 objects, have to convert to a string using np.datetime_as_string in order to set the map frame bounds.

The decorator kwargs_to_strings convert the list [w, e, s, n] to a string w/e/s/n. We need to convert datetime64 object to string internally, but it's not ideal to do the conversion in the decorator. I don't have a good solution now.

• Having even just one NaT value will raise an error, as mentioned in #242 (comment), should be a warning ideally, but it's an upstream issue. Can be tested by commenting out the 'mask' code above.

Just opened an issue (GenericMappingTools/gmt#3414) in the upstream repository.

• Support timedelta64, e.g. GPS time from a certain epoch. More of a nice-to-have really, could be done in a separate PR once this datetime one is merged in.

Perhaps we need to convert timedelta64 to double values before passing to GMT. Leave it for a separate PR.

[weiji14]

• ISO datetimes down to nanoseconds are supported by GMT!

Perhaps not true. GMT converts ISO datetime to double values internally. 2019-01-28 07:51:16.582496785 is stored as a double value 1548661876.582496881484985 (printed with %30.15lf). So GMT only supports up to microsecond, as limited by the precision of double type. Maybe I'm wrong here.

Ideally there would be a warning, but I'm not sure if it matters. The data will be plotted in pretty much the same position anyway (unless someone looks at changes every microsecond?) so it's low priority.

Just opened an issue (GenericMappingTools/gmt#3414) in the upstream repository.

Ok, I've tested the 'new' GMT master with GenericMappingTools/gmt#3415 merged in, and it's giving a warning (as expected), but the plot shows up now 🙌

Code (without mask!):

import os

os.environ["GMT_LIBRARY_PATH"] = os.path.join(os.environ["HOME"], "gmt-install-dir/lib")

import numpy as np
import pygmt

pygmt.show_versions()

# %%
x = np.array(
    [
        np.datetime64(dt)
        for dt in [
            "NaT",
            "2019-01-28T07:51:16.582496785",
            "2019-04-29T03:30:57.620034382",
            "2019-07-28T23:10:36.096639410",
            "2019-10-27T18:50:32.615134843",
            "2020-01-26T14:30:18.835180975",
        ]
    ]
)
y = np.array(
    [np.NaN, 3126.60298909, 3126.67885045, 3126.67276984, 3126.7217197, 3126.68279204,]
)

# %%
dates = np.datetime_as_string(x, unit="ns")
region = [dates[1], dates[-1], np.nanmin(y) - 0.1, np.nanmax(y) + 0.1]

# %%
fig = pygmt.Figure()
fig.basemap(
    projection="X15c/5c",
    region=region,
    frame=["WSne", "xaf+lDateTime", "yaf+lHeight(m)"],
)
fig.plot(x=x, y=y, style="t1c", pen="1p")
fig.savefig("height_over_time.png")
fig.show()

[seisman]

Passing in xarray.DataArrays to plot works (thanks to NEP18's array functionality).

@weiji14 Could you help provide a simple example for this?

And what's the best place to add some examples for calendar plots?

[weiji14]

This seems about ready, and it would be nice to merge this before we start working on #520 properly to avoid the conflicts. I'll help to fix the tests in a few minutes, seems to be just whitespace and a wrong error captured.

[weiji14]

Quick question, is there a GMT dtype for string types? I can't seem to find one at https://docs.generic-mapping-tools.org/6.1/api.html#gmt-c-api.

[seisman]

There is a GMT_TEXT, but it may not related to the GMT_Put_Strings function.

[seisman]

I'll help to fix the tests in a few minutes, seems to be just whitespace and a wrong error captured.

Thanks for help improving this PR. There are still something to do (e.g., add examples, pass a datetime region as a list), but we can merge this PR and open issues for the missing feature.

[weiji14]

There are still something to do (e.g., add examples, pass a datetime region as a list), but we can merge this PR and open issues for the missing feature.

Yep, we can add those tutorials and enhancements later 🚀