2

description of the problem

I'm writting a function to simplify my use of rasterstat.

def zonal_stat(src_rst, mask_vector, out_vector, measurments=available_stats[:4], verbose=False):
    """Compute statistics over value of the source raster for each geometry defined in the mask. results are displayed in a vector file. The user can choose the statistics he want to inclued in the final file
    
    Args : 
        src_rst (str) : path to the source raster
        mask_vector (str) : path to the vector mask 
        out_vector (str) : path to the output vector 
        measurment ([str], optional) : list of the measurment you want to perform. available measurments are :
            - min
            - max
            - mean
            - count
            - sum
            - std
            - median
            - majority
            - minority
            - unique
            - range
            - nodata
            - percentile
        the first 4 are the defaulted values. 'all' key word can also be used to perform all measurments. ['hist'] will perform an histogram of the categorical value of the raster. it will be the only stat performed if it is found in the list
        verbose (bool) : wether or not to display text
    """
    
    # apply the verbose option
    v_print = custom_print(verbose)
    
    # get the file parameters
    with rio.open(src_rst) as f:
        crs = f.crs
        nodata = f.nodata
    
    # open the mask and project it into the raster crs 
    gdf = gpd.read_file(mask_vector).to_crs(crs)
    
    categorical = False
    stats = None
    # filter the measurments
    
    if measurments: 
        if 'hist' in measurments:
            categorical = True
        elif 'all' in measurments:
            stats = available_stats[:len(available_stats) - 2 - 1]
        else:
            tmp = [stat for stat in measurments if stat in available_stats]
            stats = None if len(tmp)==0 else tmp
            
    # if categorical I need to fetch the potential values 
    if categorical:
        with rio.open(src_rst) as f:
            
            count = np.bincount(f.read(1).flatten())
            features = np.where(count!=0)[0]
            
            # remove the nodata value from the feature
            features = features[~np.isin(features, nodata)]
            
            del count
            
    
    # loop over the different geometry 
    # longer than using just rasterstats but then I cannot display any evolution informations
    tmp_res = {i:[] for i in features} if categorical else {i: [] for i in stats}

    with tqdm(total=len(gdf), disable= not verbose) as pbar:
        for index, row in gdf.iterrows():
            
            #update tqdm 
            pbar.update(1)
            
            # get the geometry coordinates 
            left, bottom, right, top = row.geometry.bounds
            
            with rio.open(src_rst) as f:
                
                # window inside the extends of the raster
                left = max(left, f.bounds.left)
                top = min(top, f.bounds.top)
                right = min(right, f.bounds.right)
                bottom = max(bottom, f.bounds.bottom)
                
                # creat the window over the geometry 
                tlx, tly = f.index(left, top)
                brx, bry = f.index(right, bottom)
                win = Window.from_slices((tlx, brx), (tly, bry))
                win_transform = f.window_transform(win)
                
                # read the data in the window
                data = f.read(1, window=win)
                
                # do a zonal stat on 1 geometry
                zs = rstats.zonal_stats(
                    [row.geometry], # there is only one geometry 
                    data, 
                    affine = win_transform, 
                    nodata = nodata, 
                    all_touched = True,
                    stats = stats,
                    categorical = categorical)[0]
                
                # release memory
                del data
                 
            # for categorical measurments I need to transform the index in int 
            if categorical:
                zs = {int(i): int(zs[i]) for i in zs}
                
            # fill the dictionary
            for idx in tmp_res:
                val = zs[idx] if idx in zs.keys() else 0
                tmp_res[idx].append(val)
        
    df = pd.DataFrame(tmp_res)
                    
    # create the result geodataframe
    stat_gdf = gpd.GeoDataFrame(
        df, 
        geometry = gdf.geometry,
        crs = crs
    )
                
    # save it 
    stat_gdf.to_file(out_vector)
    
    return df, stat_gdf

I tested this function on 2 small files and I get the following error when I run df, stat_gdf = sgt.zonal_stat(source, shapes, results, ['hist'])

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-9-c5098434ebe0> in <module>
      1 import sgt
      2 
----> 3 df, stat_gdf = sgt.zonal_stat(source, shapes, results, ['hist'])
      4 #df, stat_gdf = sgt.zonal_stat(source, shapes, results)
      5 df

~/.local/lib/python3.6/site-packages/sgt/sgt_zonal.py in zonal_stat(src_rst, mask_vector, out_vector, measurments, verbose)
    158 
    159     # save it
--> 160     stat_gdf.to_file(out_vector)
    161 
    162     return df, stat_gdf

/usr/local/lib/python3.6/dist-packages/geopandas/geodataframe.py in to_file(self, filename, driver, schema, index, **kwargs)
    744         from geopandas.io.file import _to_file
    745 
--> 746         _to_file(self, filename, driver, schema, index, **kwargs)
    747 
    748     def set_crs(self, crs=None, epsg=None, inplace=False, allow_override=False):

/usr/local/lib/python3.6/dist-packages/geopandas/io/file.py in _to_file(df, filename, driver, schema, index, mode, crs, **kwargs)
    253             crs_wkt = crs.to_wkt("WKT1_GDAL")
    254         with fiona.open(
--> 255             filename, mode=mode, driver=driver, crs_wkt=crs_wkt, schema=schema, **kwargs
    256         ) as colxn:
    257             colxn.writerecords(df.iterfeatures())

/usr/local/lib/python3.6/dist-packages/fiona/env.py in wrapper(*args, **kwargs)
    398     def wrapper(*args, **kwargs):
    399         if local._env:
--> 400             return f(*args, **kwargs)
    401         else:
    402             if isinstance(args[0], str):

/usr/local/lib/python3.6/dist-packages/fiona/__init__.py in open(fp, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt, **kwargs)
    272             c = Collection(path, mode, crs=crs, driver=driver, schema=this_schema,
    273                            encoding=encoding, layer=layer, enabled_drivers=enabled_drivers, crs_wkt=crs_wkt,
--> 274                            **kwargs)
    275         else:
    276             raise ValueError(

/usr/local/lib/python3.6/dist-packages/fiona/collection.py in __init__(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, ignore_fields, ignore_geometry, **kwargs)
    163             elif self.mode in ('a', 'w'):
    164                 self.session = WritingSession()
--> 165                 self.session.start(self, **kwargs)
    166         except IOError:
    167             self.session = None

fiona/ogrext.pyx in fiona.ogrext.WritingSession.start()

AttributeError: 'int' object has no attribute 'encode'

Can someone explain me what this error means and how to write a workaround ?

what I've already tested

I've tried on my side to create the geodataframe from the df var outside of the function (in a notebook) and there were no problem.

truc = gpd.GeoDataFrame(df, geometry=geo.geometry, crs=crs)

I've also tried to modify the code in my function with :

# create the result geodataframe
stat_gdf = gpd.GeoDataFrame(
    df, 
    geometry = gdf.geometry,
    columns = [str(i) for i in tmp_res] + ['geometry'], # geopandas refuse int as column key
    crs = crs
)

But then all the value in the geoDataFrame are replaced by Nan (None in the saved .shp)

to reproduce my sources files

def fake_tif():
    
    # initialte parameters 
    crs = rio.crs.CRS.from_epsg(4326)
    res = 0.00026949458523585647 # ~30m in deg
    bands, rows, cols = shape = (1, 6, 6)
    west, south, east, north = -cols*res/2, -rows*res/2, cols*res/2, rows*res/2
    transform = rio.transform.from_bounds(west, south, east, north, cols, rows)

    # create the dataset 
    data = np.zeros(shape, dtype=np.uint8)
    data[0, 2:4, 2:4] = 10
    data[0, 4, 4] = 10
    data[0, :2, :3] = 22
    data[0, 0, 5] = 36
    
    plt.imshow(data[0], extent=[west,east,south,north])
    plt.show()
    
    # burn it into a file 
    file = tmp_dir.joinpath('source.tif')
        
    kwargs = {
        'driver': 'GTiff', 
        'dtype': 'uint8', 
        'width': shape[1], 
        'height': shape[2], 
        'count': bands, 
        'crs': crs, 
        'tiled': False, 
        'compress': 'lzw', 
        'interleave': 'band',
        'transform': transform,
        'nodata': 0
    }
        
    with rio.open(file, 'w', **kwargs) as dst:
        dst.write(data)
        
    return file

source = fake_tif()
def fake_shape():
    """create a shahpaefile composed of 4 squares, that cross the tif file"""
    
    
    crs = rio.crs.CRS.from_epsg(4326)
    res = 0.00026949458523585647 # ~30m in deg
    
    cols = 6
    rows = 6
    west, south, east, north = -cols*res/2, -rows*res/2, cols*res/2, rows*res/2
    height = cols*res
    length = cols*res
    
    data = {
        'fixed_id' : [i+1 for i in range(4)],
        'geometry' : [
            sg.Point(west + length/4,south + 3*height/4),
            sg.Point(west + 3*length/4,south + 3*height/4),
            sg.Point(west + length/4,south + height/4),
            sg.Point(west + 3*length/4,south + height/4)
        ]
    }
    gdf = gpd.GeoDataFrame(data, crs='EPSG:4326')
    gdf.geometry = gdf.geometry.buffer(cols*res/4, cap_style=3)
    
    dst = tmp_dir.joinpath('shapes.shp')
    
    gdf.to_file(dst)
    
    gdf.plot(cmap='inferno')
    
    return dst

shapes = fake_shape()
1
  • Good reproducible example! An improvement would be to explicitly specify the imports. import rasterio as rio; import numpy as np; from matplotlib import pyplot as plt; import shapely.geometry as sg; import geopandas as gpd; from tqdm import tqdm; from rasterio.windows import Window; import rasterstats as rstats; import pandas as pd; – Stefan Jan 22 at 11:19
3

Reading this SO answer and fiona/ogrext.pyx it seems that you are passing wrong parameters to the gdf.to_file() method. According to ogrext.pyx, e.g. path, driver and crs are encoded using utf-8.

In your specific case, the problem is that your GeoDataFrame has ints as column names:

Index([10, 22, 36, 'geometry'], dtype='object')

fiona doesn't like that. A solution would be to cast the column names to str:

tmp_res = {str(i):[] for i in features} if categorical else {i: [] for i in stats}

The column names would be like this then:

Index(['10', '22', '36', 'geometry'], dtype='object')

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.