[GIS] Cannot upscale raster with rasterio

pythonrasterrasterioresolution

I am trying to upscale a raster using rasterio but I don't know why my process gets killed.

It can print "i am arrived here" but on the reproject() I have as terminal output:

killed

This is the code I am trying to run:

with rasterio.open(mask) as src:
    df_mask = src.read(1)
    aff_mask = src.transform
    profile_mask = src.profile
    width = src.width


profile_esm.update( compress='lzw', dtype = 'float32' )

## using reproject with the same CRS but changing the cell size

new_df= np.empty( shape = ( round(df_mask.shape[0] * 10),
                round(df_mask.shape[1] * 10 )), dtype = 'float32')

newaff = Affine(aff_mask.a / 10, aff_mask.b, aff_mask.c,
                aff_mask.d, aff_mask.e / 10, aff_mask.f)


profile_mask.update( compress='lzw', dtype = 'float32')

print("i am arrived here")
reproject(
        df_mask, new_df,
        src_transform = aff_mask,
        dst_transform = newaff,
        src_crs = src.crs,
        dst_crs = src.crs,
        resampling = Resampling.nearest)

print("i am after")
with rasterio.open('resampling_mask.tif', 'w', **profile_mask) as dst:
    reproject(
        df_mask, new_df,
        src_transform = aff_mask,
        dst_transform = newaff,
        src_crs = src.crs,
        dst_crs = src.crs,
        resampling = Resampling.nearest)

    dst.write(new_df,  indexes=1)

Best Answer

I'm not sure why you're seeing the death of your program, but the reproject function is not the way to upscale data. Instead, use the dataset's read method and specify a larger output shape:

with rasterio.open(mask) as dataset:
    upscaled_mask = dataset.read(1, out_shape=(dataset.height * 10, dataset.width * 10))

Related Solutions

Rasterio – Reading Raster Files by Block

Reading Rasters by block can be done in rasterio and I'd argue it's easier than in GDAL. There is even a tutorial on windowed read/write over at GitHub.

Let's take a look at the read functions arguments, which allows you to set a window to read data from:

def read(self, indexes=None, out=None, window=None, masked=False,
        boundless=False):
"""Read raster bands as a multidimensional array
    Parameters
    ----------
    window : a pair (tuple) of pairs of ints, optional
        The optional `window` argument is a 2 item tuple. The first
        item is a tuple containing the indexes of the rows at which
        the window starts and stops and the second is a tuple
        containing the indexes of the columns at which the window
        starts and stops. For example, ((0, 2), (0, 2)) defines
        a 2x2 window at the upper left of the raster dataset.
"""

So all you need now is a way to create the window tuples. That's what rasterios block_windows function is for:

def block_windows(self, bidx=0):
    """Returns an iterator over a band's block windows and their
    indexes.
    [....]
    The primary use of this function is to obtain windows to pass to
    read_band() for highly efficient access to raster block data.
    """

So in your case I guess the code would look something like this:

with rasterio.open("your/data/geo.tif") as src:
     for block_index, window in src.block_windows(1):
         block_array = src.read(window=window)
         result_block = some_calculation(block_array)

[GIS] Mask rasterio raster with geopandas shapefile

If your geodataframe has a column named "geometry" and geometry consists of Shapely Polygons, the following works for me

from rasterio.mask import mask

out, _ = mask(data, gdf.geometry, invert=False)

Related Question

python – Updating Dimensions of Raster with Rasterio Without Scaling Data