Working with point data

[adamjstewart]

A user asked me how to use point data in TorchGeo, so I'm writing a summary of how I would do this. Note that this is experimental and subject to change.

I think the workflow would look something like this:

# Define datasets
sentinel = Sentinel2(root="...")
gbif = GBIF(root="...")
dataset = sentinel & gbif

# Define sampler and data loader
sampler = PointGeoSampler(image=sentinel, point=gbif, size=256, length=10000)
dataloader = DataLoader(dataset, batch_size=128, sampler=sampler, collate_fn=stack_samples)

# Iterate over mini-batches, train models like you would in torchvision
for batch in dataloader:
    image = batch["image"]
    bbox = batch["bbox"]

The workflow is almost identical to https://github.com/microsoft/torchgeo#geospatial-datasets-and-samplers but with one key difference: PointGeoSampler.

PointGeoSampler hasn't been written yet, but it would behave something like this. First, instead of taking in a single dataset as input, it would take in image and point datasets separately. The __init__ method would look something like:

hits = []
for p in point:
    if p in image:
        hits.append(p)
self.hits = hits

Basically, we're just trying to make sure all of our points are within the bounds of the images in our image dataset.

The __iter__ method would look something like:

for _ in range(len(self)):
    # pick a random point
    # pick a random bounding box around that point
    yield bounding_box

I haven't yet found the time to work on PointGeoSampler. There's a research project I'm working on that may need this eventually, so I'm hoping to implement this later this fall. But if anyone wants to take a stab at this, I would be happy to review a PR!

[ak3ra]

Hi @adamjstewart,

I am looking into this now and would like to start thinking about the implementation of the PointGeoSampler.

Say I have an instance of the GBIF dataset, which contains lat/lon points would the PointGeoSampler take this as well as say Sentinel2Dataset separately?

or are you assuming the PointGeoSampler already should contain the satellite image and the lat/lon points?

Let me know your thoughts

[adamjstewart]

Yes, GBIF and Sentinel2 would be separate datasets and PointGeoSampler would take in 2 datasets as input instead of 1.

[ak3ra]

Hi @adamjstewart, could you please take a look at this code and tell me if this is in the right direction for this? or what is missing?

The PointGeoSampler takes in two datasets, a point dataset e.g (GBIF, INaturalist, EDD), and a Raster Dataset e.g Sentinel2.

I would like to understand how:

1. The Sampler would return an intersectional patch containing a point (lat/lon) given if the implementation of GBIF, INaturalist datasets index bounding boxes, not 2d lat/lon coordinates.

class PointGeoSampler(GeoSampler):
  """Samples from an intersection of points and image GeoDatasets.

  This sampler should take in a pair of coordinates and region of interests and
  returns image patches containing the coordinates of interest.
  """

  def __init__(
      self,
      dataset: GeoDataset,  ## for exampler raster dataset of sentinel 2
      points: GeoDataset,   ## e.g GBIF dataset containing coordinate points 
      roi: Optional[BoundingBox] = None,
      shuffle: bool = False,
  )-> None:
    super().__init__(dataset, roi)
    self.shuffle = shuffle
    self.points = points
    self.dataset = dataset

    self.hits = []

    for p in self.points.index.bounds:
      if p in dataset.index.bounds:
        self.hits.append(p)

  def __iter__(self) -> Iterator[BoundingBox]:
    generator: Callable[[int], Iterable[int]] = range
    if self.shuffle: 
      generator = torch.randperm 

    for idx in generator(len(self)):
      yield BoundingBox(*self.hits[idx].bounds)

    def __len__(self) -> int: 
      return len(self.hits)

[ak3ra]

I just made a PR Please review let me know if this is along the lines of what you suggest.

[adamjstewart]

Does it work?

[ak3ra]

It is a work in progress, one issue is that the two datasets do not have a matching crs. so the intersection of the two Rtrees do not really overlap, trying to fix that.

[adamjstewart]

That one is actually easy to fix:

points.crs = dataset.crs

[ak3ra]

That fixes it.
Now I need to test it with an overlapping dataset. I am trying to query Gbif data within the Chesapeake region so I can test it with the tutorial NAIP dataset