Splitting raster data into equal pieces with GDAL in Python

[@louis-philip]

Thanks for all this, it helps a lot!

[@t0mmy-Sun]

quality content!! massive thanks

[@albertnkwasa2648]

This channel is just dope! Thank you👌

I tested the code and I found two issues:
1) using gdal warp i got a Warning: "Warning 1: for band 1, destination nodata value has been clamped to 0, the original value being out of range." It seems related with the nodata values. At the end the resulting tiles are not workables.
2) the resolution res is not valid for y-axis, therefore the tiles generated are incomplete. I defined resx and resy (one per axis), at the result was ok.
But in general, your code is the best example that I found surfing the web. Thanks a lot.

[@joegawlik9901]

You're saving my life with these tutorials🥲❤

[@jasonbushell7080]

Your series of videos has really helped me, so thank you! I have split a 100MB file into 9, and each piece is now over 1GB! I'm not sure what I've done wrong? 🤣

[@louis-philip]

Same here! Individual clipped images are about 8 times the file size of the original full res image.

[@ryanmurphyri]

SO GOOD!

[@johnowusukonduah2305]

Any idea about how to do deep learning with geotiff files?

[@jaiswalfelipe1269]

Thank you! Solved my problem. However, does this work on multispectral images?

[@jaiswalfelipe1269]

I tried to follow your method on 4000x4000x23 tif file, I used div = 16, I was expecting to output 256 patches with 250x250x23 in shape (16,000,000 / 62500 = 256 patches), but It gave me 250 patches instead. What could be the reason for this? the reso? Thank you in advance!

[@kerimmirzeyev4340]

Hello, thank you for this tutorial it helped me a lot of, I just have a 1 question, I modfied your code instead of 1 data, I'm trying to divide folder of sentinel 2 images to 32x32 tiles, even though I don't get error and when printing it does show the coordinates of newly divided tiles, but output folder remains empty I think warp is not doing what it supposed to do and I don't know what else I could to do get it working, can you tell me what I'm doing wrong?
This is my code
import os
from osgeo import gdal
# Specify the input folder containing 256x256 images
input_folder = r"E:\ytu\PHD\data
eslihan_sentinel_testing\kodlar\SmartpolAIProcess\images"
# Output folder for 32x32 tiles
output_folder = r"E:\ytu\PHD\data
eslihan_sentinel_testing\imagesTiles"
image_files = [file for file in os.listdir(input_folder) if file.endswith('.tif')]
# Loop through each image file
for image_file in image_files:
sen = gdal.Open(os.path.join(input_folder, image_file))

# Get the geotransform
gt = sen.GetGeoTransform()
# Get the dimensions of the image
xlen = gt[1] * sen.RasterXSize
ylen = gt[1] * sen.RasterYSize
# Number of tiles in x and y direction
div = 8
# Size of a single tile
xsize = xlen // div
ysize = ylen // div
# Create lists of x and y coordinates
xsteps = [gt[0] + xsize * i for i in range(div + 1)]
ysteps = [gt[3] - ysize * i for i in range(div + 1)]
# Loop over min and max x and y coordinates
for i in range(div):
for j in range(div):
xmin = xsteps[i]
xmax = xsteps[i + 1]
ymax = ysteps[j]
ymin = ysteps[j + 1]

print("xmin: "+str(xmin))
print("xmax: "+str(xmax))
print("ymin: "+str(ymin))
print("ymax: "+str(ymax))
print("
")

# Use gdal warp to create 32x32 tiles
gdal.Warp(
os.path.join(output_folder, f"{os.path.splitext(image_file)[0]}_{i}_{j}.tif"),
sen,
outputBounds=(xmin, ymin, xmax, ymax),
dstNodata=-9999,
)

[@JuanCarlosMH]

Hi! Thanks for this video & your channel!! I have a photogrammetry produced raster with local coordinates, dem.GetProjection() displays 'LOCAL_CS["Local Coordinates (m)",UNIT["metre",1,AUTHORITY["EPSG","9001"], when I run the gdal.Warp command It results in ERROR 1: PROJ: proj_uom_get_info_from_database: unit of measure not found. Any ideas??

[@JuanCarlosMH]

Somehow worked with gdal.Translate method

[@gwynrivers1219]

Hi, how/where do I specify a folder for the newly generated DEMs?

[@jeppert7215]

Maybe this comes too late, but: You can add the folder path in front of the string in the gdal.Warp or gdal.Translate functions. So for example, you could write gdal.Warp("output/dem" + str(i) + str(j)... etc). I hope it still helps :)

[@gwynrivers1219]

@@jeppert7215 Yes, I figured this out, but thank you for the reply!

[@KamranAli-kx6rb]

I am converting rasters into 5 * 5 patches using the gdal. The script is working fine but some patches which are at the boundary their dimensions that are not correct. Dimensions are 5*5 but there are 2,3 pixels not 5 * 5 dimension pixels other pixels are no data pixels. I want to split patches if they show 5 * 5 dimensions there should be 5 * 5 dimension pixels. patches at boundary if there is no proper 5 * 5 dimension patch it should show exact dimensions it should not show 5 * 5 dimensions.
Code:
import os, gdal
in_path = 'D:\\im\\'
input_filename = 'image.tif'
out_path = 'D:\\split2\\'
output_filename = 'tile_'
tile_size_x = 5
tile_size_y = 5
ds = gdal.Open(in_path + input_filename)
band = ds.GetRasterBand(1)
xsize = band.XSize
ysize = band.YSize
for i in range(0, xsize, tile_size_x):
for j in range(0, ysize, tile_size_y):
com_string = "gdal_translate -of GTIFF -srcwin " + str(i)+ ", " + str(j) + ", " +
str(tile_size_x) + ", " + str(tile_size_y) + " " + str(in_path) +
str(input_filename) + " " + str(out_path) + str(output_filename) + str(i) + "_" +
str(j) + ".tif"
os.system(com_string)