Reading QA_Pixel Channel in Landsat8

[TolgaAktas]

As part of my dataset instantiation, QA_PIXEL is also a band that I am reading along with other 11 channels in L8 products. QA_PIXEL is supposed to be a number of bit masks but it's currently being read as float values, which I am not sure if the Dataset object is reading correctly. How can I correctly read the QA_PIXEL band as part of band stack?

bands = ["B1","B2","B3","B4","B5","B6","B7","B8","B9","B10", "B11","QA_PIXEL"] cloudfree_dataset = Landsat8(root=cloudfree_path,bands=bands)

[adamjstewart]

#1149 may help with this, although we'll likely need to go a step further and not cast anything ever.

[TolgaAktas]

They shouldn't have any bounding box errors really, right?

[adamjstewart]

Bounding box errors == no intersection

That shouldn't happen as long as you modified filename_regex to exclude the date and there are indeed images with spatial overlap in both datasets.

[TolgaAktas]

they are literally from the same site. I am thinking that this might be the problem. We created different versions of datasets, whereas this one is the one where we cropped the inner part of the dataset (I believe via GDAL?) to avoid sampling empty pixels. With the original images datasets, the intersection does work, but we sample nonpixels too often. Could you tell me what we should be paying attention in cropping the images, are we missing a key information in the cropped images that goes into the bbox calculation?

Regardless, I think they should still be intersecting, since this is literally the same area. But as I said, it's either using the modified regex and getting identical bands, or not using it and thus not having a successful intersection.

[adamjstewart]

I ran your script with the data you sent me using the latest version of main. The script completed without errors, and shadow_img and clean_img are not identical. Maybe a different set of images is required to reproduce the issue? Or maybe the issue is in your plotting code?

[adamjstewart]

Also, answers to your questions over email:

What would happen if I tried UnionDataset instead of Intersection dataset?

You would get one larger dataset that contains both cloudy and cloud-free images. Sampling this dataset would randomly choose one or the other, you wouldn't get both at the same time.

For intersecting some datasets, where there’s no separation of clean and shadow image datasets, the size of the intersection dataset is the minimum of the sizes of the two datasets. When I modify the regex as you suggested, the size of the intersection dataset is the multiplication of the sizes of the two datasets.

The length of an IntersectionDataset is the number of possible overlaps between images in the two datasets. For example, if dataset 1 has 10 files, and dataset 2 has 20 files, and there is only 1 file in each dataset that overlaps, then an IntersectionDataset would have length 1. The size of that image would be the area of intersection between the two overlapping images.

[TolgaAktas]

As it stands, the sample code I sent you couldn't find the bands so I had to make the following edition:

1. Overriding bands & rgb_bands for the Landsat8 class.
2. Removed the bands argument from the Landsat8 constructor.

root_path = "/dirs/data/euphrates/datasets/landsat/summer_trimmed_copy"
bands = ["B1","B2","B3","B4","B5","B6","B7","B8","B9","B10", "B11","QA_PIXEL"]
rgb_bands = ['B4', 'B3', 'B2']

Landsat8.all_bands = bands
Landsat8.rgb_bands = rgb_bands
 
shadow_image_path = os.path.join(root_path,"shadow_image")
clean_image_path = os.path.join(root_path,"clean")
shadow_image = Landsat8(root=shadow_image_path)
clean_image  = Landsat8(root=clean_image_path)

At the bottom, I added the following checks to see if they are identical bands:

Landsat8.filename_regex = ".*_(?P<band>[A-Z0-9_]+)\."
shadow_image = Landsat8(root=shadow_image_path,bands=bands)
clean_image  = Landsat8(root=clean_image_path,bands=bands)
all_set = shadow_image & clean_image
sampler = RandomBatchGeoSampler(all_set,size=(512,512),length=100,batch_size=10)
dataloader = DataLoader(all_set,batch_sampler=sampler,collate_fn=stack_samples)
for elem in dataloader:
    a = elem
    break
batch = a["image"]
shadow_img = batch[0,:12,:,:]
clean_img  = batch[0,12:,:,:]
print(torch.allclose(shadow_img[0],shadow_img[1]))
print(torch.allclose(clean_img[0],clean_img[1]))
print(torch.allclose(shadow_img,clean_img))

Allclose within separate images return True, although the last allclose returns False. So all the bands within a single image is still identical.

1. I think it doesn't have anything to do with the plotting function since allclose already shows that they are identical bands
2. You mentioned that it works without any errors on the latest version of main. I installed torchgeo via "pip install torchgeo==0.4.0", do you think that might be the reason ? What fixes have gone into the latest main that might not be in 0.4?

[TolgaAktas]

Welcome back home :) Feel free to drop by RIT sometime! I wanted to try a minimally viable solution for the home town first, plus we get a good variety of clouds here. What do you think might be the problem? I am going through geo.py and landsat.py to pinpoint. Couple interesting observations:

1. self.filename_glob = self.filename_glob.format(bands[0]) in landsat.py. I bet when it's reading identical bands, it's always the first band.
2. I think the problem might lie in the merge_dataset-indices where it's trying to find the hits. I can successfully print a number of hit1 and hit2 bboxes, but I don't know much about rtree lib so I couldn't investigate further.
3. The error message pops up when I try to print the dataset object, which prompts str function. It should print self.bounds, but I couldn't find anywhere in the code where self.bounds gets set to a value. I see that the minx maxx are printed equal which shouldn't happen, and indeed I can print out a number of bounds for the datasets. But the scale of the minx and maxx are unrealistically high (10 to the 308), which tells me that they both get assigned to a maximum value somewhere as a default?

[adamjstewart]

Fixed: #1191

[adamjstewart]

Thanks for finding and reporting all these bugs! This one has been in TorchGeo since before our first release!!

[TolgaAktas]

Thank you so much for your continued support and patience with it! This is super exciting since we have been trying to resolve it for a while, and now this will facilitate some really interesting work, but also that I am glad to contribute to such a big open-source project for the first time as torchgeo! Hope I can help more!

That being said, I think it might be still worthy to look into why the IntersectionDataset fails to find the intersecting indices when we don't alter the regex string, since the datasets are from the same scene, and they don't get matched up. Any ideas as to why that might be happening ? (especially for the minx and maxx pair, where they get set to an extremely large number ~ e+308? )

Again, thank you very much :)

[adamjstewart]

There's two comments I would make about that.

First, while these images overlap spatially, they don't overlap temporally, so the intersection does not exist. All datasets in TorchGeo keep track of both the spatial and temporal location of each image. So that's why there's no intersection unless you tell the base class to ignore the time component.

Second, the very large numbers you are seeing are the largest and smallest floating point numbers in python. When you try to compute the bounds of an R-tree containing no entries, rtree returns these bounds. There's a bug in rtree that results in these bounds being invalid. If you want to try to fix that bug, I can point you to the issue. But for your case, all it means is that IntersectionDataset is empty because there is no spatiotemporal overlap.

Hope that helps clarify the behavior you're seeing.