I'm trying to reproduce the R rasterize function in python to streamline my code and am having issues with summing multiple points in the raster. Essentially I want to sum the value of my variable for every point in the grid I create. However currently it only appears to be returning the value of the first point in the grid it sees.

import rasterio
import geopandas as gpd

def create_raster_defined_res(geodataset, output, out_var, resolution):  
    res = (resolution, resolution)
    bounds = geodataset.total_bounds
    kwargs = {
            'count': 1,
            'crs': geodataset.crs['init'],
            'bounds': bounds,
            'width': max(int(ceil((bounds[2] - bounds[0]) / float(res[0]))), 1),
            'height': max(int(ceil((bounds[3] - bounds[1]) / float(res[1]))), 1),
            'affine': Affine(res[0], 0, bounds[0], 0, -res[1], bounds[3]),
            'transform': (bounds[0], res[0], 0, bounds[3], 0, -res[1]),
            'compress': 'lzw',
            'driver': u'GTiff',
            'dtype': 'float64',
            #'nodata': -1e+308

    with rasterio.open(output, 'w', **kwargs) as out:
        out_arr = np.empty([kwargs['height'],kwargs['width']])

        # this is where we create a generator of geom, value pairs to use in rasterizing
        shapes = ((geom,value) for geom, value in zip(geodataset['geometry'], geodataset[out_var]))

        burned = rasterize(shapes=shapes, fill=0, out=out_arr, transform=kwargs['affine'])
        out.write_band(1, burned)

create_raster_defined_res(geopandas_dataset, output_tiff, output_var, 2500)

Is there a way to add a parameter to sum the value of all the points that enter into the grid or tell it to add instead of just taking the value of the first point?

  • Hej! Did you find a solution for your problem?? I'm having kind of the same here. :) – Christin Abel Feb 25 at 17:21
  • Yup I eventually found out that the default function for rasterize is to replace and not add. Look at the Merge_alg property in rasterize -> rasterio.readthedocs.io/en/stable/api/… – DarknessFalls Feb 26 at 19:18
  • merge_alg=MergeAlg.add – DarknessFalls Feb 26 at 19:19
  • Oh and I I created my own indicator = 1 per record in a prior step which I then added. My application required % so I ran a raster layer for each category of my variable and then summed the fields to calculate the %. EX: n_50 = src_50.read(1, window = Window(0,i*100, width, 100)) low_pct_50 = np.true_divide(src_50.read(2, window = Window(0,i*100, width, 100)),n_50) – DarknessFalls Feb 26 at 19:21

Your Answer

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

Browse other questions tagged or ask your own question.