Topics

rasterio.windows.transform seems to not scale my windows correctly, am I using it wrong?


Ryan Avery
 

I have a landsat band that I have windowed so that I now have about 200 512x512 windows in a list called chip_list_full. I have plotted these on top of my aoi to confirm that the windowing worked and found that it did not work as expected following this stack overflow example code where I computed a transform for each window using the original dataset transform. the result is that windows are spatially, diagonally offset from the upper left corner of the original image



each window in the image above was plotted like so, with a unique transform called custom_window_transform:
chips_with_labels = []
fig, ax = plt.subplots()
gdf.plot(ax=ax)
for chip in chip_list_full:
    custom_window_transform = windows.transform(chip[0], band.transform)
    window_bbox = coords.BoundingBox(*windows.bounds(chip[0], custom_window_transform))
    window_poly = rio_bbox_to_polygon(window_bbox)
    gpd.GeoDataFrame(geometry=gpd.GeoSeries(window_poly), crs=band.meta['crs'].to_dict()).plot(ax = ax, cmap='cubehelix')

When I change the custom_window_transform to instead be the original band transform, the result is correct



chips_with_labels = []
fig, ax = plt.subplots()
gdf.plot(ax=ax)
for chip in chip_list_full:
    window_bbox = coords.BoundingBox(*windows.bounds(chip[0], band.transform))
    window_poly = rio_bbox_to_polygon(window_bbox)
    gpd.GeoDataFrame(geometry=gpd.GeoSeries(window_poly), crs=band.meta['crs'].to_dict()).plot(ax = ax, cmap='cubehelix')

where "band.transform" is my original image transform, and the dataset I have windowed.

My question is, what is the purpose of windows.transform and is there some other way I should be using it in this case? Or is my use of the original dataset transform correct?

Thanks for the feedback!
where "band" is 


Sean Gillies
 

Hi Ryan,

I've been on vacation, just now getting the time to answer questions. Answers below.

On Wed, Aug 21, 2019 at 3:13 PM Ryan Avery <ravery@...> wrote:
I have a landsat band that I have windowed so that I now have about 200 512x512 windows in a list called chip_list_full. I have plotted these on top of my aoi to confirm that the windowing worked and found that it did not work as expected following this stack overflow example code where I computed a transform for each window using the original dataset transform. the result is that windows are spatially, diagonally offset from the upper left corner of the original image



each window in the image above was plotted like so, with a unique transform called custom_window_transform:
chips_with_labels = []
fig, ax = plt.subplots()
gdf.plot(ax=ax)
for chip in chip_list_full:
    custom_window_transform = windows.transform(chip[0], band.transform)
    window_bbox = coords.BoundingBox(*windows.bounds(chip[0], custom_window_transform))
    window_poly = rio_bbox_to_polygon(window_bbox)
    gpd.GeoDataFrame(geometry=gpd.GeoSeries(window_poly), crs=band.meta['crs'].to_dict()).plot(ax = ax, cmap='cubehelix')

When I change the custom_window_transform to instead be the original band transform, the result is correct



chips_with_labels = []
fig, ax = plt.subplots()
gdf.plot(ax=ax)
for chip in chip_list_full:
    window_bbox = coords.BoundingBox(*windows.bounds(chip[0], band.transform))
    window_poly = rio_bbox_to_polygon(window_bbox)
    gpd.GeoDataFrame(geometry=gpd.GeoSeries(window_poly), crs=band.meta['crs'].to_dict()).plot(ax = ax, cmap='cubehelix')

where "band.transform" is my original image transform, and the dataset I have windowed.

My question is, what is the purpose of windows.transform and is there some other way I should be using it in this case? Or is my use of the original dataset transform correct?

Thanks for the feedback!

Your use of windows.transform in the first case is correct, and your use of windows.bounds in the second case is correct. Where you went wrong, and probably due to sketchy documentation, is in

    window_bbox = coords.BoundingBox(*windows.bounds(chip[0], custom_window_transform))

The second argument of the windows.bounds function (see https://rasterio.readthedocs.io/en/latest/api/rasterio.windows.html#rasterio.windows.bounds) must be the affine transformation for the "dataset" on which we're applying the given window, not any other affine transformation. I really must explain this better in the docs, sorry about that.

--
Sean Gillies