Re: Create raster aligned with one created by gdal_translate

thomaswanderer@...
 
Edited

I had a very similar approach to calculate the new affine & width/height manually as you. But I have now replaced it with just rasterio functions and when comparing the results of the new affine, they were identical :-)

I just wonder, why you add 0.5 when calculating the height & width? My code looked like this (windows is the extent I am matching my raster to, profile is the target raster profile, copied from my source raster and then modified)
if window is not None and window.is_polygon:
                    # Make sure we operate with the same coordinates
                    window = window.transform(profile['crs'].to_epsg())

                    # Get the coordinates from raster and window bounds
                    s_xmin, s_ymin, s_xmax, s_ymax = raster_bounds.bounds
                    w_xmin, w_ymin, w_xmax, w_ymax = window.bounds

                    # Get the source pixel width and height
                    x_pixel = profile['transform'].a
                    y_pixel = profile['transform'].e

                    # Calculate the shift vector from source bounds to target bounds
                    shift_xmin = w_xmin - s_xmin
                    shift_ymin = w_ymin - s_ymin
                    shift_xmax = w_xmax - s_xmax
                    shift_ymax = w_ymax - s_ymax

                    # Normalize the shift vector by the source pixel width and height
                    shift_xmin_norm = ceil(shift_xmin / x_pixel) * x_pixel
                    shift_ymin_norm = floor(shift_ymin / y_pixel) * y_pixel
                    shift_xmax_norm = floor(shift_xmax / x_pixel) * x_pixel
                    shift_ymax_norm = ceil(shift_ymax / y_pixel) * y_pixel

                    # Create the new boundaries of the target raster
                    new_xmin = s_xmin + shift_xmin_norm
                    new_ymin = s_ymin + shift_ymin_norm
                    new_xmax = s_xmax + shift_xmax_norm
                    new_ymax = s_ymax + shift_ymax_norm

                    # Calculate width & height (Pixel) for the target raster
                    new_width = abs(int((new_xmax - new_xmin) / x_pixel))
                    new_height = abs(int((new_ymax - new_ymin) / y_pixel))
                    if new_width == 0 or new_height == 0:
                        raise Exception(f'New raster has invalid measure (width or height is 0 pixels)')

                    # Update profile with windows dimensions
                    profile.update(width=new_width)
                    profile.update(height=new_height)

                    # Create new affine and update profile with it
                    new_affine = rasterio.transform.from_origin(
                        new_xmin,
                        new_ymax,
                        new_width,
                        new_height
                    )
                    profile.update(transform=new_affine)

But i have replaced all this now by:


# Make sure we operate with the same coordinates
# Make sure we operate with the same coordinates window = window.transform(profile['crs'].to_epsg()) # Create a new window and affine from the window extent raster_window = source.window(*window.bounds) raster_window = raster_window.round_offsets(op='ceil').round_lengths(op='floor') # Update target raster with the new affine new_affine = source.window_transform(raster_window) profile.update(transform=new_affine) # Update target raster with new windows shape (Make sure raster stays within extent) c_xmin, c_ymin, c_xmax, c_ymax = source.window_bounds(raster_window) w_xmin, w_ymin, w_xmax, w_ymax = window.bounds x_size, y_size = source.res corrected_width = floor(abs(c_xmin - w_xmax) / x_size) corrected_height = floor(abs(w_ymin - c_ymax) / y_size) profile.update(width=corrected_width) profile.update(height=corrected_height)

Both methods actually created and affine with the same values, that is why I replaced my own code with the rasterio functions.

Join main@rasterio.groups.io to automatically receive all group messages.