Resample raster from one resolution to another with GDAL Python API

Question

i have some trouble warping one raster to different resolution. My input is a one channel tiff with. I know the geo-coordinates of x_min and y_max. I also know that the pixel resolution is 5 meters. I now want to warp the same raster to pixel resolution of 6.5 meters.

Input in the function is the numpy array and the upper left corner coordinates with the src and dst pixel resolution.

My code is:

def resample_raster(raster, raster_x0, raster_ymax, src_xres, src_yres, dst_xres, dst_yres):
 dst_height=int((dst_yres/src_yres)*raster.shape[0])
 dst_width=int((dst_xres/src_xres)*raster.shape[1])
 dstAR=np.empty((dst_height, dst_width))
 srcDS=dstDS=gdal_array.OpenArray(raster)
 srcDS.SetGeoTransform((raster_x0,src_xres,0, raster_ymax, 0, -src_yres ) )
 dstDS=gdal_array.OpenArray(dstAR)
 srcDS.SetGeoTransform((raster_x0,dst_xres,0, raster_ymax, 0, -dst_yres ) )
 gdal_warp_options=osgeo.gdal.WarpOptions(xRes=dst_xres,yRes=dst_yres, resampleAlg=osgeo.gdal.GRA_NearestNeighbour)
 osgeo.gdal.Warp(dstDS, srcDS,options=gdal_warp_options)
 new_raster=np.array(dstDS.GetRasterBand(1).ReadAsArray())
 return new_raster

The function gives me an empty array or is killed.

Answer 1

I am not certain of what you try to achieve but I have two sugestions:

1. Follow this idea on how to resize a numpy array using scipy.
2. Use gdal.Warp to warp a tiff to a different resolution. For example:

import gdal
infn = '/path/to/source.tif'
outfn = '/path/to/target.tif'
xres=6.5
yres=6.5
resample_alg = 'near'
ds = gdal.Warp(outfn, infn, warpoptions=dict(xRes=xres, yRes=yres, resampleAlg=resample_alg)
ds = None

Answer 2

I did this with the command line gdalwarp from my python script as follows:

# Build command 
command = f'gdalwarp -t_srs EPSG:{targetEPSG} -r cubic -tr {targetResolution} {targetResolution} -te {leftBound} {bottomBound} {rightBound} {topBound} "{inputRasterPath}" "{newPath}"'
 
# Reproject vector 
os.system(command)

Answer 3

I have downsampled a worldview-3 image (GeoTIFF format) to a WorldView-2 resolution, using gdal/osr Linux package.

First, I implement a special line to take the resolution of a reference image (image with lower resolution):

X_RES = gdalinfo -mm $PATH_TO_REFERENCE_IMAGE_WITH_LOWER_RESOLUTION | grep "Pixel Size" | grep "\([0-9]\.[0-9]\+\),円\([-w][0-9]\.[0-9]\+\)"

Than, applied the gdalwarp routine to each image in a input folder. The complete solution:

PATH_TO_HIGHER_RESOLUTION_IMAGES=1ドル
PATH_TO_REFERENCE_IMAGE_WITH_LOWER_RESOLUTION=2ドル
PATH_TO_OUTPUT=3ドル
RESAMPLE_METHOD="cubic"
X_RES="0.5"
# TODO: get pixel size
#X_RES = gdalinfo -mm $PATH_TO_REFERENCE_IMAGE_WITH_LOWER_RESOLUTION | grep "Pixel Size" | grep "\([0-9]\.[0-9]\+\),円\([-w][0-9]\.[0-9]\+\)"
#echo $X_RES
if [ -z "$X_RES" ]; then
 echo "\$X_RES resolution is empty"
fi
for image in $PATH_TO_HIGHER_RESOLUTION_IMAGES*.TIF; do
 filename=$(basename $image)
 if [[ ! -e $PATH_TO_OUTPUT ]]; then
 mkdir $PATH_TO_OUTPUT
 elif [[ ! -d $PATH_TO_OUTPUT ]]; then
 echo "$PATH_TO_OUTPUT already exists but is not a directory" 1>&2
 fi
 result_image=$PATH_TO_OUTPUT$filename
 gdalwarp -of GTiff -wo NUM_THREADS=ALL_CPUS -co TILED=YES -co TILED=YES -co COMPRESS=LZW -tr $X_RES $X_RES -r $RESAMPLE_METHOD $image $result_image
done

and usage:

./resample.sh ../folder-with-higher-resolution/ ../folder-with-lower-resolution/REFERENCE_IMAGE.TIF ../output-resampled-images/

** Observe the -tr and -ts flags. First one is for pixel size, where another is for desired image dimension.