Time Series Support - Test and Use Cases

[nilsleh]

This discussion is a continuation of #877 but could also serve for further use cases so that we get a better understanding of the different use cases people might desire for Time Series tasks.

[nilsleh]

Hi @Spiruel, could you check how many hits you are getting for the geolocations here maybe that can help to see whether the sampling issue is in the spatial or the time domain. This is also a good reminder that I should include error messages or warnings that explain why one would see the error you are seeing.

Also looking at the code snippet you posted, it looks like you are passing only the input_ds and now the ds to the sampler. Is it possible that the input and target dataset have different locations or projections? The ForecastDataset inherets from IntersectionDataset and ensures there is overlap between the datasets.

[Spiruel]

Hi @nilsleh thank you for your help looking into this.

I debugged the hits and saw there was an issue in the time domain. Turns out I didn't have the correct time_unit argument for the month of images I was looking in. I'm not getting the previous error anymore (I'm getting spatial and temporal hits), but I am coming across other issues eg RTreeError: Coordinates must be in the form (minx, miny, maxx, maxy) or (x, y) for 2D indexes.

I left in the ForecastDataset in my code snippet but I don't expect to use it for my use case of classifying/segmenting time series. What I'm looking for is the SequentialGeoSampler to yield available input time series of shape (num_observations, channels, width, height) with an associated target label. Could you let me know what the correct configuration should be to achieve this, and I can work from there to debug, offer feedback?

I've attached some example input data examples.zip - there should be chips with multiple temporal hits in that month.

Thanks

#prechipped harmonised landsat-sentinel rasters  across a calendar year
class TimeSeriesInputRaster(RasterDataset):
    filename_glob = "*.tif"
    filename_regex = r"^\d+_\d+_(?P<date>\d{8})T?.tif$"
    date_format = "%Y%m%d"
    is_image = True
    separate_files = False
    all_bands = ["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B09", "B10", "B11", "B12"]

#yearly cover crop masks for oct-apr
class TimeSeriesTargetRaster(RasterDataset):
    filename_glob = "*_cc2class.tif"
    filename_regex = r"""
        ^(?P<date>\d+)
        _cc2class.*$
    """
    date_format = "%Y"
    is_image = True
    separate_files = False
    all_bands = ["target"]

input_root = "S30/examples"
label_root = "data"

input_ds = TimeSeriesInputRaster(input_root, as_time_series=True)
label_ds = TimeSeriesTargetRaster(label_root, as_time_series=True)

ds =  ForecastDataset( #not using this as not a forecast problem
    input_dataset=input_ds,
    target_dataset=label_ds
)

sampler = SequentialGeoSampler(
    input_ds,
    size=32,
    length=3,
    encoder_length=1,
    prediction_length=1,
    time_unit="hours", #this has to be set to hours to work. Shouldnt it be days?
)

sample = input_ds[next(iter(sampler))]
# input = sample["image"] # [sample_window, num_bands, height, width]
# target = sample["target"] # [target_window, num_targets, height, width]

#expect input features of size [3, 14, 32, 32] and corresponding label of size [1, 2, 32, 32]

[Spiruel]

Hi @nilsleh unfortunately I can't share my targets data and I don't have it prechipped... But if you use cdl = CDL(root="cdl_data", download=False, checksum=True) on a maize or soy class then they should work absolutely fine for a MWE.

Thanks for the clarification about SequentialGeoSampler, that makes sense.

[Spiruel]

By num_observations I mean the time dimension. So for example in other semantic segmentation tasks in agriculture eg. torchgeo.datasets.ZueriCrop yields a fixed number 71 observations per sample (shape is 71, 9, 24, 24 I think).

[nilsleh]

Another question, which is a very general question that we have not decided on how to handle it is the fact that most geospatial time-series data is irregularly spaced. The tif files given in your zip folder also illustrate this:

                       filename       date
0   501900_4476600_20191008.tif 2019-10-08
1   502860_4455480_20191008.tif 2019-10-08
2   502860_4455480_20191013.tif 2019-10-13
3   501900_4476600_20191013.tif 2019-10-13
4   501900_4476600_20191015.tif 2019-10-15
5   561420_4484280_20191015.tif 2019-10-15
6   502860_4455480_20191015.tif 2019-10-15
7   501900_4476600_20191023.tif 2019-10-23
8   502860_4455480_20191023.tif 2019-10-23
9   502860_4455480_20191124.tif 2019-11-24
10  561420_4484280_20191124.tif 2019-11-24
11  501900_4476600_20191124.tif 2019-11-24
12  501900_4476600_20191204.tif 2019-12-04
13  502860_4455480_20191204.tif 2019-12-04
14  561420_4484280_20191204.tif 2019-12-04
15  502860_4455480_20191207.tif 2019-12-07
16  501900_4476600_20191207.tif 2019-12-07
17  501900_4476600_20200101.tif 2020-01-01
18  502860_4455480_20200101.tif 2020-01-01
19  502860_4455480_20200202.tif 2020-02-02
20  561420_4484280_20200202.tif 2020-02-02
21  501900_4476600_20200202.tif 2020-02-02
22  501900_4476600_20200420.tif 2020-04-20
23  502860_4455480_20200420.tif 2020-04-20
24  501900_4476600_20200427.tif 2020-04-27
25  561420_4484280_20200427.tif 2020-04-27
26  502860_4455480_20200427.tif 2020-04-27

For example if the encoder length is specified to be 4 with a time unit of "days", and you get a sample that looks for hits between 2019-10-13 and 2019-10-17, there will only be 2 hits and the returned sample will have a time dimension of 2. This problem is what I am continuing to work on for the general implementation in the PR. However, maybe I have also misunderstood your usecase and we can find a solution for that.

[Spiruel]

In response to your comment about irregularly spaced time series - I suppose for my use case I'm interested in all available hits to be returned, and random subsampling being used to ensure a fixed number of observations in the final sample.

[nilsleh]

Okay, so if I understand correctly if you have a query that corresponds to a Sampler with encoder_length=4 and time_unit=days then your expectation for 2 example queries is the following:

1. from 2019-10-13 to 2019-10-17 return 2 time dimension hits that are merged geospatially, so you get [2, 14, 32, 32] ([num_time_hits, channels, x, y] as your sample dimension and a target that corresponds to the latest available time point to form a target [1, 2, 32, 32] ([time, num_targets, x, y])
2. from 2020-04-20 to 2020-04-24 returns 1 time dimension hit that is merged geospatially to form [1, 14, 32, 32] as the sample dimension and a target of [1, 2, 32, 32]

If this is what you want, then you will need to define another collate_fn for a dataloader that combines the individual samples to a batch because by default pytorch expects each sample to have the same dimensions, which cannot be guaranteed.

Or could you elaborate how you expect the random sub sampling procedure to guarantee a fixed number of observations if the time-series is irregularly spaced and at odds with the current implementation of SequentialGeoSampler that is looking for a fixed time-horizon only.

Just to give another example, if for that same data shown in the dataframe, the encoder_length=2 and time_unit=month (so asking for data in a 2 month time horizon) the following could happen:

1. from 2019-10-08 to 2019-12-07 so 2 month period there are 7 time units so you would get a tensor of dimension [7, 14, 32, 32]

I agree that the naming scheme is not ideal at the moment to show that the encoding_length is connected to the time horizon one specifies and not a fixed tensor length being returned.

[Spiruel]

If you have encoder_length=4 and time_unit=days then I follow the two example queries you give. A target corresponding the latest time point would be great and would be what I require for my use case. If we can get a MWE of that working then I can run with that to code up a few time series classification examples.

I think having this option, with a collate_fn to handle the subsampling is a good design choice. Like this approach in Breizhcrops. It is what I assumed SequentialGeoSampler did before I played around with the code.

Yes the name encoding_length did confuse me to begin with, but I think I get it now!

Thank you :)

[nilsleh]

Thanks for the example of Breizhcrops. I am not familiar with that dataset but from briefly looking at the code it also appears that the dataset is pre_chipped so it is previously defined what a single sample is and that it consists of a fixed time length and also the spatial dimension because the dataset is a pytorch map-style Dataset that just retrieves samples.

The use case that I am currently working on in the PR has the assumption that you have a collection of tif files covering a certain area and possibly multiple files for the same region. However, samples are not predefined beforehand so through a sampler we want to generate random samples both in the spatial and time domain. The current SequentialGeoSampler is intended to generate samples where you have a input sequence of a certain length and the target sequence follows consecutively. This is intended to be used for a forecasting task from historical data towards the future. However, there are of course many different time-series prediction tasks that are different (even if only slightly) in their sampling logic and we hope to cover these different tasks with different GeoSamplers.

I only started with the SequentialGeoSampler since that is a use case I have for a project, but your use case sounds interesting and quiet general so I will look into how we can get your example to work.

[Spiruel]

Thank you - if you have time to look at this then that would be appreciated. I have to admit that I initially thought SequentialGeoSampler was quite generalised and I could use it for my use case out of the box.

However if this is not the case, then getting this working would be well worth it - my use case is quite ubiquitous and would be useful to many if it's integrated into torchgeo.

[nilsleh]

Hi @Spiruel I have been working on a new Sampler design that is a bit more universal. There are still no unit tests for it and still subject to change, however, with the committed updates to that branch, the following example works for me:

from torchgeo.datasets import RasterDataset, MultiQueryDataset, CDL

from torchgeo.samplers import TimeWindowGeoSampler

input_root = "your_example_data/examples"
label_root = "cdl"

class TimeSeriesInputRaster(RasterDataset):
    filename_glob = "*.tif"
    filename_regex = r"^\d+_\d+_(?P<date>\d{8})T?.tif$"
    date_format = "%Y%m%d"
    is_image = True
    separate_files = False
    all_bands = ["B01", "B02", "B03", "B04", "B05", "B06", "B07", "B08", "B8A", "B09", "B10", "B11", "B12"]

input_ds = TimeSeriesInputRaster(input_root, as_time_series=True)
label_ds = CDL(label_root, years=[2019])

ds =  MultiQueryDataset(
    input_dataset=input_ds,
    target_dataset=label_ds
)

sampler = TimeWindowGeoSampler(
    ds,
    size=32,
    length=100,
    encoder_length=3,
    prediction_length=1,
    consecutive=False,
    time_unit="months"
)

for query in sampler:
    sample = ds[query]

Besides the naming update, the new sampler argument consecutive=False has the effect that the target sequence is not afterwards in time but the last time-step in the input sequence if prediction_length=1, so I think this is what you are looking for in your use case.

I am working on different use cases now to try to find edge cases etc. where the current implementation would break and additionally try to test it with unit tests so this is still just a draft I am afraid, but maybe you can see if this is closer to what you are looking for. If you want to use this with a DataLoader you would still have to define a collate_fn because the samples generated are not guaranteed to have the same time length so automatic batching will fail. We are working on different approaches to this as well.

Let me know how this goes for your setup, happy to help and improve our implementation and explanations like naming and docstring.

[adamjstewart]

Some interesting time-series formulations from climate models (will append to this list as I find more):

• https://arxiv.org/abs/2202.11214 (section 2.2): predict next step during pre-training, predict next 2 steps during fine-tuning
• https://arxiv.org/abs/2211.02556 (section 3.4): train different models for predicting 1-hr, 3-hr, 6-hr, and 24-hr lead times, then aggregate
• https://arxiv.org/abs/2212.12794 (section GraphCast): train autoregressive models on N next steps (N varies from 1..12)

[calebrob6]

Dropping some use cases:

As a user, I'd like to be able to download all the Sentinel 2 data for a year for a single MGRS tile, then create a Sentinel2 dataset object and sample one scene at random per month. This should give me (12 months x num channels x height x width) "imagery" objects.

Further, as a user, I'd like to be able to take the above Sentinel2 dataset object, intersect it with a CDL dataset and get an annual "mask" key too (so that I can predict CDL labels from a monthly time series of S2).

In the first case, instead of "select one scene at random per month", I'd also like to be able to "median composite", "select the least cloudy scene", etc. (select least cloudy would probably need a STAC dataset).