0

I'm using rasterio to read in a .tif file and converting it to vectors using the rasterio.features module. For a particularly large .tif file, I thought I'd use the windowed read option, but for some reason the projections are coming way off when doing the windowed read. Notice how the Y coordinates are very different between the two outputs but the X coordinates are similar. Any ideas on why this is happening? See the example below:

With windowed read

with rasterio.open('raster2vector/clipped_fire2_reprojected.tif') as dataset: for block_index, window in dataset.block_windows(1): mask = dataset.read_masks(1, window=window) image = dataset.read(1, window=window) results = ( {'properties': {'pixel_value': v}, 'geometry': s} for i, (s, v) in enumerate( shapes(image=image, mask=mask, transform=dataset.transform))) geoms = list(results) if not geoms: continue print(dataset.crs) print(geoms[0]) break

produces

CRS({'init': 'epsg:2163'}) {'geometry': {'coordinates': [[(-1861517.9955081018, 916182.1706994411), (-1861517.9955081018, 915914.427799729), (-1860447.0239092538, 915914.427799729), (-1860447.0239092538, 916182.1706994411), (-1861517.9955081018, 916182.1706994411)]], 'type': 'Polygon'}, 'properties': {'pixel_value': 6.0}}

Without windowed read

with rasterio.open('raster2vector/clipped_fire2_reprojected.tif') as dataset: mask = dataset.read_masks(1) image = dataset.read(1) results = ( {'properties': {'pixel_value': v}, 'geometry': s} for i, (s, v) in enumerate( shapes(image=image, mask=mask, transform=dataset.transform))) print(dataset.crs) geoms = list(results) print(geoms[0])

produces

CRS({'init': 'epsg:2163'}) {'geometry': {'coordinates': [[(-1861517.9955081018, -457874.39062262117), (-1861517.9955081018, -458409.87642204505), (-1861250.2526083898, -458409.87642204505), (-1861250.2526083898, -458142.1335223331), (-1860982.509708678, -458142.1335223331), (-1860982.509708678, -458409.87642204505), (-1860714.7668089657, -458409.87642204505), (-1860447.0239092538, -458409.87642204505), (-1860447.0239092538, -457874.39062262117), (-1861517.9955081018, -457874.39062262117)]], 'type': 'Polygon'}, 'properties': {'pixel_value': 6.0}}

1

Turns out the transform object needs to actually be transform = rasterio.windows.transform(window, dataset.transform) for windowed reads.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.