0

I'm iterating through thousands of polygons. Each polygon clips a tif raster producing a numpy array. Then the numpy array values are manipulated and the tif raster is updated with those new values. The script works for a few thousand polygons but at some point it throws an error.

import os, sys, datetime, time
import geopandas as gpd
import rasterio
from rasterio.mask import mask
from fiona.crs import from_epsg
import numpy as np
import pycrs
from rasterio.windows import Window
from shutil import copy
import petl as etl
import shapely.wkt

alr_path = os.path.join(working_folder, 'INPUT', 'NI_COMBINED_ALR.tif')
alr = rasterio.open(alr_path)

out_meta = alr.meta.copy()   

alr_avg_path = os.path.join(working_folder, 'OUTPUT', 'alr_avg.tif')    
copy(alr_path, alr_avg_path)

feature_gdf = gpd.read_file(gdb_polygs, driver='FileGDB', layer=fc_vml_polygs)

for index, row in feature_gdf.iterrows():
    window, out_transform = mask(alr, row, all_touched=True, crop=True)

    x = alr.meta['transform'][2]
    y = alr.meta['transform'][5]

    width = window.shape[2]
    height = window.shape[1]
    col_off = out_transform[2]
    row_off = out_transform[5]

    window_orig = alr.read(1, window=Window(int((col_off - x)/5), int((y - row_off)/5), width, height))

    if np.all([window < 0]):
        avg = -1
    else:
        avg = window[window != -1].mean().item()
    window_avg = (np.where(window!=-1, int(round(avg,0)), window_orig))

    with rasterio.open(alr_avg_path, "r+", **out_meta) as dest:
        dest.write(window_avg, window=Window(int((col_off - x)/5), int((y - row_off)/5), width, height))

This is the error that the last line of the code above throws:

ERROR 1: LZWDecode:Corrupted LZW table at scanline 10808
ERROR 1: TIFFReadEncodedStrip() failed.
ERROR 1: .\OUTPUT\alr_avg.tif, band 1: IReadBlock failed at X offset 0, Y offset 10808: TIFFReadEncodedStrip() failed.
Traceback (most recent call last):
  File "footprint.py", line 277, in <module>
    process_raster(con, pg_table_name, path_tif)
  File "footprint.py", line 206, in process_raster
    dest.write(numpy_array, window=Window(x, y, width, height))
  File "rasterio\_io.pyx", line 1409, in rasterio._io.DatasetWriterBase.write
OSError: Read or write failed

Any ideas about what this error means and how it could be overcome?

1 Answer 1

-1

The error hints at malformed or broken raster input data. rasterio is unable to decompress some data. Try if you can open the raster file for which this error occurs in QGIS.

Your code seems to be fine.

3
  • Polygon feature class can be opened with no problems in ArcGIS. Polygon object where error is thrown is read successfully, it clips the raster producing an array successfully and only when updating raster file with this array the error is thrown. Jul 24, 2019 at 12:29
  • I meant the raster file, sorry for not being clear earlier. Please add a minimal complete code example that triggers the problem to your question text. Jul 24, 2019 at 13:24
  • Ah ok, you meant the raster. I can't see any problem when opening the raster file in ArcGIS. I've added a minimal complete code to the description. Jul 24, 2019 at 14:38

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.