Reproject, resample and clip raster data with GDAL in Python

In this tutorial, I explain how to use gdalwarp in Python to reproject raster data to a different coordinate reference system, change the resolution (resample) and clip it to a shapefile.
GDAL/OGR Python API: gdal.org/python/
Code:
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
ds = gdal.Open("dem.tif")
reproject
dsReprj = gdal.Warp("demReprj.tif", ds, dstSRS = "EPSG:4326")
resample
dsRes = gdal.Warp("demRes.tif", ds, xRes = 150, yRes = 150,
resampleAlg = "bilinear")
clip
make sure your raster data and shapefile have the same projection!
dsClip = gdal.Warp("demClip.tif", ds, cutlineDSName = "star.shp",
cropToCutline = True, dstNodata = np.nan)
visualize
array = dsClip.GetRasterBand(1).ReadAsArray()
plt.figure()
plt.imshow(array)
plt.colorbar()
plt.show()
close your datasets!
ds = dsClip = dsRes = dsReprj = None

[@908572813]

you are simply umbelievable! great job i would like to spend hours and hours listening and seeing your videos. many tnx

[@drforest]

Your videos have been so helpful to me. Thank you. Your videos are helping us track salmon habitat and monitor regional urbanization. 🎉

[@sayyadsahilahmed321]

These are really helpful and to the point, I got all what I need here. Many thanks, keep educating us👍

[@Antonio-lp8hx]

Awesome, I absolutely love your channel :)

[@qasim4296]

@1:10 Reproject @3:48 Resample @5:37 Clip I’m afraid you will stop making these valuable videos once you are done with your PhD. Please keep making, if possible

[@CarlosSánchez-m9m]

Great tutorial!! I'm stuck on this, and I was wondering if you could help me out. I have a coarse mesh coarse.tif and a fine mesh fine.tif. The extent of coarse.tif is larger and completely includes fine.tif. Both have different resolutions, but I want the points (pixels) of each to be "nested", giving precedence to the arrangement of points in the fine mesh. This means that starting from the fine mesh, I want to construct a coarse mesh where there is one coarse mesh point for every 4 fine mesh points. In other words, the square formed by every 4 fine mesh points has the corresponding coarse mesh point in the center. Once this coarse mesh is created, the average of the involved fine mesh points is considered in the points where both meshes overlap, and interpolation of the coarse mesh points is used outside the overlap. Is it possible by using gdalwarp? Thanks

[@yusranurtanrverdi7101]

This tutorial saved my day. Thank you so much

[@malequinn]

Great tutorial, I am however still looking for how to make the plot be on the coordinate plane

[@geetarNL]

Awesome, thanks :) do you know any way of feeding gdal.Warp a multipolygon that is a gdal.geometry object to clip the input image? I used CreateGeometryFromWkt to create a gdal geometry multipolygon, but Warp does not accept this for the cutline.

[@anantharamaniyer9135]

Happy New Year to you. Please spare some time in perhaps writing a blog post for us to refer to your wonderful lessons. How much of a difference in file size does it make if we go from a 25m resolution to a 100m resolution - would it be a factor of 4?

[@makingsenseremotely6207]

I don’t have time for blogposts, sorry. I am currently super busy with my PhD project. When I downsample my dem.tif file from 25m to 100m, its size is reduced by factor 15.

[@XinhLe]

Hi, Thanks so much for the tutorial. Please make a video on how to use different methods of georeferencing/registration to georeference or register an image to a set of control points?

[@martinsz441]

very nice video, thanks!

[@engr.abdullahazzamsafi4515]

Hi, Really you are amazing. It's a great job. You didn't mentioned the script how to save it after resampling. Moreover, How we can resample a multiple rater images at once. I will be grateful to comment me the script. Thanking you in anticipation

[@makingsenseremotely6207]

Hi, you don't need to save the resampled raster - gdal.Warp will automatically create a new file (demRes.tif) - check your working directory. To resample multiple images, simply list all the files that you want to resample and process them in a loop.

[@engr.abdullahazzamsafi4515]

@@makingsenseremotely6207 Thank you so much. Can you please provide any script for calculating RMSE of multiple raster files and obtain spatially RMSE map. I will be very thankful to you. Thank you

[@vijaygarg13]

Hi, Can we use color map in gdal. I have a NDVi image but it is in greyscale.Can we change the image from greyscale to single band rgb image?

[@DianaRocíoGalindoGonzález]

Graciaaaaas

[@TiriAlain]

Thank you for your tuto ! Can you help me ? I must to clip an area around a coordonate (GPS or another) from a hudge .tif (in fact 2, a DSM and a DTM). This area can be a circle or a polygone of 20 or 25 meters. It's seem to be similar like you are doing here.

[@melissalopezportillopurata1134]

Thank you so much! I always find the solution to my questions on your videos.

[@brettkaiser4467]

sorry to be offtopic but does anyone know of a way to log back into an Instagram account?? I was stupid forgot my login password. I would love any assistance you can give me!

[@colinbruce3171]

@Brett Kaiser Instablaster :)

[@brettkaiser4467]

@Colin Bruce thanks so much for your reply. I found the site on google and Im trying it out now. Takes quite some time so I will reply here later with my results.

[@brettkaiser4467]

@Colin Bruce it did the trick and I finally got access to my account again. Im so happy! Thanks so much, you saved my ass !

[@colinbruce3171]

@Brett Kaiser Glad I could help xD

[@AlessioDiLorenzo]

Great tutorial!

[@mdkashif2560]

Hi, many thanks for sharing this valuable video! Could you please guide me as to how I can clip with gdal from a raster image when I’ve got only four coordinates? I don’t have a shape file. 🙏

[@nesrinekadri2597]

Thanks for you. I want to know how to clip band 4 and 8 of many sentinel 2 imagery ( shapefile area)

[@RAJATGARG45]

Thanks for a great tutorial. I was trying to resample my raster data but I got following error: ERROR 1: Attempt to create 0x0 dataset is illegal, sizes must be larger than zero. My Code: src1 = gdal.Open("E:/fire/Similipal/Code/Data/Pre.tif") src1_ = gdal.Warp("XYZ.tif", src1, xRes = 20, yRes =20) Any idea how to resolve this? Thanks.

[@makingsenseremotely6207]

My guess is that your input raster uses a geographic CRS in units of degrees. The output file resolution must be specified in target georeferenced units. As you are not attempting to change the CRS, gdal will assume you want an output resolution of 20x20° which is probably much larger than your dataset, thus the error. So either figure out what 20m is in degrees or change your CRS to a projected one that records locations in meters.

[@anushriyajain1959]

Hii, thanks for the video. :) I am working with different datasets from different satellites and reanalysis products to be used in single model. So when i am using gdal to reproj, resamp, and clippping all datasets to same extent, i am getting 1 or 2 extra rows in few datasets. So when they are kept above each other their spatial extent is coming little different. Could you help me how to tackle this problem? i tried with geopandas as well but same issue. Thanks in advance.

[@makingsenseremotely6207]

Could be related to different resolutions of your datasets, especially if you clip them first and resample afterwards. If, for example, you want your clipped dataset to have an extent of xmin=0 and xmax=100, that will work for a 10m image, but if your dataset has a x-resolution of 3m, you won’t be able to exactly fit all pixels in there and will get a clipped image with a max. x-extent of either 99 or 102. So make sure that your data all have the same resolution before you start clipping.

[@anushriyajain1959]

@@makingsenseremotely6207 Thank you so much. I am doing that only as you suggested but still same issue. Still Thanks alot

[@anantharamaniyer9135]

Hello again - I tried resampling from 25 to 100m resolution and got the error message - AttributeError: type object 'object' has no attribute '__getattr__' Any ideas how to get around this?

[@makingsenseremotely6207]

This error message usually appears when you try to get an attribute that doesn’t exist, e.g. when you run ds.SomeNonExistingAttribute. Maybe you misspelled something, like GetRasterBand?

[@FocusStudySoundscapes]

Data Set available?