resampling a geotiff to smaller cell size

Sam McClatchie

unread,

Jul 7, 2023, 1:52:58 AM7/7/23

to WhiteboxTools

Hello

I am using Whitebox Workflows for python on Ubuntu Linux. I want to resample a geotiff of sea surface temperature with approximately 1 km resolution to a smaller cell size and save the resampled raster. The input raster was geolocated with gdal translate and gdal warp.

Here is the gdal info on the input raster:

-----------------------------------------------------------------

$ gdalinfo Kaipara_Manukau_NZ_overlay_proj.tif

Driver: GTiff/GeoTIFF
Files: Kaipara_Manukau_NZ_overlay_proj.tif
Size is 1400, 872
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["unknown"],
AREA["World - 85°S to 85°N"],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (19247139.958156999200583,-4328176.075298279523849)
Pixel Size = (143.148641046410006,-143.148641046410006)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
TIFFTAG_XRESOLUTION=300
TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left (19247139.958,-4328176.075) (172d54' 0.00"E, 36d12' 0.00"S)
Lower Left (19247139.958,-4453001.690) (172d54' 0.00"E, 37d 5'58.63"S)
Upper Right (19447548.056,-4328176.075) (174d42' 1.07"E, 36d12' 0.00"S)
Lower Right (19447548.056,-4453001.690) (174d42' 1.07"E, 37d 5'58.63"S)
Center (19347344.007,-4390588.883) (173d48' 0.53"E, 36d39' 4.05"S)
Band 1 Block=1400x1 Type=Byte, ColorInterp=Red
Mask Flags: PER_DATASET ALPHA
Band 2 Block=1400x1 Type=Byte, ColorInterp=Green
Mask Flags: PER_DATASET ALPHA
Band 3 Block=1400x1 Type=Byte, ColorInterp=Blue
Mask Flags: PER_DATASET ALPHA
Band 4 Block=1400x1 Type=Byte, ColorInterp=Alpha

---------------------------------

This is my Python script:

import os
import rasterio
import matplotlib.pyplot as plt
import imageio

import whitebox
wbt = whitebox.WhiteboxTools()

# Pass in the region list variable from shell script
from sys import argv
region_str = argv[1]
region_list = region_str.split(",")


for region in region_list:

    input_raster = imageio.v2.imread("/mnt/data/dynamic_data/projects/projects2022/development_recreational_fishing_high_resolution/figures_Avenza/"+region+"_overlay_proj.tif")
    # input_raster = wbt.read_raster("/mnt/data/dynamic_data/projects/projects2022/development_recreational_fishing_high_resolution/data/" + region + "_sst_reference.tif")

    # print("########### input raster #################")
    print(f'Rows: {input_raster.configs.rows}')
    print(f'Columns: {input_raster.configs.columns}')
    print(f'Resolution (x direction): {input_raster.configs.resolution_x}')
    print(f'Resolution (y direction): {input_raster.configs.resolution_y}')
    print(f'North: {input_raster.configs.north}')
    print(f'South: {input_raster.configs.south}')
    print(f'East: {input_raster.configs.east}')
    print(f'West: {input_raster.configs.west}')
    print(f'Min value: {input_raster.configs.minimum}')
    print(f'Max value: {input_raster.configs.maximum}')
    print(f'EPSG code: {input_raster.configs.epsg_code}') # 0 if not set
    print(f'Nodata value: {input_raster.configs.nodata}')

    # NOT USED: 
    # base_raster = wbt.read_raster("/mnt/data/dynamic_data/projects/projects2022/development_recreational_fishing_high_resolution/data/" + region + "_topo_reference.tif")

    # i = "/mnt/data/dynamic_data/projects/projects2022/development_recreational_fishing_high_resolution/figures_Avenza/"+region+"_overlay_proj.tif"
    out = "/mnt/data/dynamic_data/projects/projects2022/development_recreational_fishing_high_resolution/figures/SST_plot_" + region +"_resampled_raster.tif"

    resampled_raster = wbt.resample(
       [input_raster],  
       [out],
       cell_size=10,
       base=None,
       method="cc",
    )

And this is the error that is returned:

-------------------------------------------------------

>>> resampled_raster = wbt.resample(
... [input_raster],
... [out],
... cell_size=10,
... base=None,
... method="cc",
... )
./whitebox_tools --run="Resample" --inputs='[array([[[129, 129, 129, 255],
[129, 129, 129, 255],
[129, 129, 129, 255],
...,
[128, 127, 127, 255],
[130, 129, 129, 255],
[129, 129, 129, 255]],

[[129, 129, 129, 255],
[133, 131, 129, 255],
[150, 139, 128, 255],
...,
[148, 121, 106, 255],
[160, 135, 123, 255],
[136, 130, 127, 255]],

[[130, 130, 129, 255],
[137, 132, 129, 255],
[169, 147, 127, 255],
...,
[125, 76, 48, 255],
[137, 92, 69, 255],
[133, 123, 117, 255]],

...,

[[130, 129, 129, 255],
[138, 132, 130, 255],
[175, 143, 135, 255],
...,
[148, 175, 167, 255],
[143, 167, 160, 255],
[132, 138, 137, 255]],

[[130, 129, 129, 255],
[133, 130, 129, 255],
[148, 135, 132, 255],
...,
[125, 137, 133, 255],
[124, 134, 131, 255],
[119, 122, 121, 255]],

[[129, 129, 129, 255],
[129, 129, 129, 255],
[129, 129, 129, 255],
...,
[ 98, 98, 98, 255],
[ 98, 98, 98, 255],
[ 98, 99, 99, 255]]], dtype=uint8)]' --output='['/mnt/data/dynamic_data/projects/projects2022/development_recreational_fishing_high_resolution/figures/SST_plot_Kaipara_Manukau_NZ_resampled_raster.tif']' --cell_size='10' --method=cc -v --compress_rasters=False

****************************
* Welcome to Resample *
* Powered by WhiteboxTools *
* www.whiteboxgeo.com *
****************************
Reading data...
thread 'main' panicked at 'called `Option::unwrap()` on a `None` value', whitebox-raster/src/lib.rs:1361:69
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

-----------------------

So far I have not been able to run this code with the RUST BACKTRACE, sorry.

The input geotiff file can be accessed here:

http://u.pc.cd/aKNctalK

The python script can be accessed here. If you want to run it, region = "Kaipara_Manukau_NZ"

http://u.pc.cd/a31rtalK

Any help would be gratefully accepted.

Sam

Whitebox Geospatial Inc

unread,

Jul 7, 2023, 9:12:39 AM7/7/23

to Sam McClatchie, WhiteboxTools

Hello Sam,

Your script looks good but there is one major problem with it. You are attempting to input a rasterio (GDAL) in-memory raster and are expecting the wbt.resample function to return one, but this isn’t the case. The WhiteboxTools Python API is really just a thin wrapper around the WhiteboxTools command line application, which takes command line arguments as inputs and outputs. This means that all inputs and outputs must be specified as file path & name strings. Here is an example of how you would write this script using WhiteboxTools:

import whitebox

wbt = whitebox.WhiteboxTools()

wbt.set_working_dir("/path/to/data/") # Sets the working directory

wbt.resample(

inputs="my_raster.tif",

output="my_resampled_raster.tif",

cell_size=10,

method="cc"

)

print("Done!")

The description of the Resample tool contains a description of the Python API for the tool here:

https://www.whiteboxgeo.com/manual/wbt_book/available_tools/image_processing_tools.html#Resample

There is also a very important document outlining how to use Whitebox with Python in the user manual here:

https://www.whiteboxgeo.com/manual/wbt_book/python_scripting/scripting.html

If you do in fact need to deal with in-memory rasters, so that you can read the properties of the raster, specify in-memory rasters as inputs to tools, and have those tools output in-memory rasters, then what you really need is Whitebox Workflows for Python (WbW) instead of WhiteboxTools.

https://www.whiteboxgeo.com/whitebox-workflows-for-python/

Unlike WbT, WbW is not open-source and costs $10 / user. Generally, because WbW is a native Python extension library, the equivalent WbW scripts are much more natural than the equivalent WbT script. For example, the equivalent resample WbW script would look like this:

import whitebox_workflows

wbe = whitebox_workflows.WbEnvironment('Geomorphometry2023')

wbe.set_working_dir("/path/to/data/") # Sets the working directory

raster = wbe.read_raster('my_raster.tif') # create an in-memory Whitebox raster object

# Get some information about the raster

print(f'Rows: {raster.configs.rows}')

print(f'Columns: {raster.configs.columns}')

print(f'Resolution (x direction): {raster.configs.resolution_x}')

print(f'Resolution (y direction): {raster.configs.resolution_y}')

print(f'North: {raster.configs.north}')

print(f'South: {raster.configs.south}')

print(f'East: {raster.configs.east}')

print(f'West: {raster.configs.west}')

print(f'Min value: {raster.configs.minimum}')

print(f'Max value: {raster.configs.maximum}')

print(f'EPSG code: {raster.configs.epsg_code}') # 0 if not set

print(f'Nodata value: {raster.configs.nodata}')

# Run the resample function

resampled = wbe.resample(

inputs=[raster],

cell_size=10,

method="cc"

)

# Save the file to disc

wbe.write_raster(resampled, "resampled_raster.tif", compress=True)

print("Done!")

WbW offers the ability to read in rasters (as well as vectors and lidar data) and to retrieve raster properties and even to manipulate individual grid cells. Take a look at the ‘Working with raster data’ section of the WbW user manual for more information:

https://www.whiteboxgeo.com/manual/wbw-user-manual/book/working-with-raster-data.html

I hope that heps.

Regards,

John

--
You received this message because you are subscribed to the Google Groups "WhiteboxTools" group.
To unsubscribe from this group and stop receiving emails from it, send an email to whiteboxtool...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/whiteboxtools/35ed8128-b948-48eb-902c-4afbbd135f68n%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Whitebox Geospatial Inc

unread,

Jul 7, 2023, 9:26:10 AM7/7/23

to Sam McClatchie, WhiteboxTools

Hi Sam,

I see now in re-reading your email, that in fact you do have a license for WbW already. Okay, with that info in mind, then I would say that the issue is that you are using the WbT library (‘import whitebox’) rather than the Whitebox Workflows library (‘import whitebox_workflows’). And, as I said in the previous response, the two libraries result in very different types of scripts. You will need to use something that looks more like the second script in my response. Also, you don’t need to use rasterio to read in rasters and get their properties in your script, you can do it all from Whitebox Workflows instead. That’s one of the joys of using WbW. Again, for a reference and guidance on how to write WbW scripts, I’d really suggest looking over the WbW user manual, which I find quite useful. Best of luck.

Cheers,

John

From: whiteb...@googlegroups.com <whiteb...@googlegroups.com> on behalf of Sam McClatchie <smccl...@fishocean.info>
Date: Friday, July 7, 2023 at 1:53 AM
To: WhiteboxTools <whiteb...@googlegroups.com>
Subject: resampling a geotiff to smaller cell size

--

Sam McClatchie

unread,

Jul 8, 2023, 12:11:12 AM7/8/23

to WhiteboxTools

Thanks John, that was extraordinarily helpful. The adjustments to the import statements quickly fixed my problem. And thanks for the link to the Manual for WbW. I look forward to beginning to use WbW regularly in my work.