Read and write raster files with GDAL in Python
ฝัง
- เผยแพร่เมื่อ 20 ก.ย. 2024
- This tutorial explains how to read raster data as an array and save arrays as a GeoTiff file using the GDAL library in Python.
GDAL/OGR Python API: gdal.org/python/
Code:
from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
import
ds = gdal.Open("dem.tif")
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
band = ds.GetRasterBand(1)
array = band.ReadAsArray()
plt.figure()
plt.imshow(array)
manipulate
binmask = np.where((array \> = np.mean(array)),1,0)
plt.figure()
plt.imshow(binmask)
export
driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("binmask.tif", xsize = binmask.shape[1],
ysize = binmask.shape[0], bands = 1,
eType = gdal.GDT_Int16)
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(binmask)
outband.SetNoDataValue(np.nan)
outband.FlushCache()
close your datasets and bands!!!
outband = None
outds = None
I want cry of happyness, I´m so grateful with you, finally i have found what I as looking for, i has been a great tutorial.Greetings from Perú. :3
This has been very helpful, I remember taking a course of Python in GIS and being very lost however I continued. Following this tutorial allowed me to finally make sense of those codes I've been using because they worked somehow. So big thank you!
Gracias. Me sacó de un gran apuro el hecho de exportar un raster filtrando valores específocos.
Definitely appreciate the tutorial!
Thanks for making this series you saved me
Wow!!! You did a very great and accurate lesson, thank you very much, dear.
Hoping you will do more in time series raster volume spatial analysis.
Great tutorial, thanks! The third and fifth values in gt are skew
The zeros in the geotransform refer to rotation. Zero in this case because north is up
Thank you !
Great great great video. Thank you for doing these, they’re so helpful.
Thanks for sharing these great learning resources!
Muchas gracias, clarísimo el video, entendí todo!
Very useful and informative!
This is amazing, thank you!
Thanks so much.
Well explained, thank you!
Obrigado pela playlist, Otimo Conteudo !!
thank you,this video is very very very nice! give you a big like!
Awesome tutorial :)
Great channel - thank you!
thanks
It's Great. Thank you
I'm using python to compute the average of multiple rasters and save it as another raster (to use in QGIS) and the raster I end up with only contains the values 0 and 1 when the array is storing floating point values. Other than that issue, great tutorial, thank you!
I changed the eType from Int32 to Float32 and that fixed the issue
Very clear tutorial. Many thanks for making this seem so simple ! I was wondering if there are any benefits to converting the processed tiff to a CSV for further manipulation via pandas or dask etc or is there a reason for keeping the output in tiff ?
If you want to do further processing with pandas, you can definitely go ahead and turn the raster into CSV format, however, CSV files with XY-coordinates and corresponding values for every cell tend to get really big if you work with larger raster data sets. Also, the geographic information (projection, coordinate system, etc.) won’t be stored. As I am using GDAL a lot for raster data processing and manipulation, keeping my output in tiff format is certainly easier, as GeoTIFF is one of the supported raster drivers and CSV not.
Very great and clear work you have done Ma'am. i have a question please, i am working on my raster dataset for prediction like ANN, RF and CNN, i have converted the rasters into numeric and then train and test the data and got very good accuracy. Now i need to convert my test data into raster again as final prediction map but i don't know how to do this, please guide me thanks.
Hi @Channel For Everyone, Channel For Everyone, I am also working on similar project with raster dataset of rainfall data, Would you mind walking me through on how you approached the problem and how you were to convert the data to numeric and carry out ML training on it. Please reply, thanks.
I think you can write it to raster using gdal and numpy. You should know number of rows and columns and reshape it and write it. You can copy projection etc from your original raster.
Thanks for your contribution. It is really useful!
I have two questions.
How can I import a RasterStack as an array with 3 dimensions?
How can I import the RasterStack or RasterBand files into a folder as a list of arrays?
Thank youu for yout time!!!
thx!it's very helpful!!!!!
hi, thanks for those videos. are you plan to make video for gdal in python on vector geoprocessing too?
Hi, I'm definitely planning to make more videos about gdal in Python, and I'm always open to suggestions on the specific topic, so vector data is now on the list :)
RasterIO or GDAL??
what is the better to learn? since I have finished Pandas and GeoPandas course.
Thanks
Thanks you !
When I run the script to test gdal.Open, no variable shows up under variable explorer. How can I fix this please?
I want to generate the raster data, so what should be the method i should follow should be it on python or any other ?
Please suggest me some
I am using PyCharm and it doesn't display the two boxes that you have on the right half of the screen (the table and the Console 1/a) so I don't get to see what I'm executing.
Thank you!!! Very usefull videos, and I love the sound of your keyboard... Hahaha what model is?
Haha, thanks. It's just the built-in one from my mac.
thanks for this useful video , I have a question , how we can change the resolution of raster file ?
I am getting an error: Unable to allocate 3.36 GiB for an array with shape (59622, 60471) and data type bool
If I want to use no data values instead of zeros for line 19 would I need to move line 32 above line 19? outband.SetNoDataValue (np.nan)
binmask = np.where((array >= np.mean(array)),1,np.nan)
What if I want multiple bands in my output, will I have to run a loop and store all the bands in the output file. I am dealing with hyperspectral data.
I had trouble writing results from 3 bands numpy array to geotif.
One band is working. It would be great if you have example of writing from numpy array to multibeand image (RGB for example).
Thanks.
Here is my code modified to create a raster with 3 bands, each storing the binmask array. Hope that helps.
# export
driver = gdal.GetDriverByName("GTiff")
driver.Register()
outds = driver.Create("rgb.tif", xsize = binmask.shape[1],
ysize = binmask.shape[0], bands = 3,
eType = gdal.GDT_Int16)
outds.SetGeoTransform(gt)
outds.SetProjection(proj)
for i in range(1,4):
outband = outds.GetRasterBand(i)
outband.WriteArray(binmask)
outband.FlushCache()
outband = None
outds = None
What if I got a crop of the numpy array and would like to keep the global coordinates of the remaining pixels? Is there some way to pass the Affine Transform or the geo location of some of the array corners to the driver?
Thanks!
How is your IDE connected to your QGIS Software, I am curious to know how to make this happen, Does it support other IDE's as well?
Thank you very much for your tutorial. I followed all the steps, but when I open at ArcGIS, my tiff has value from 1 to 1. It should be 0 or 1, because is a binary mask like yours. Do you have any tip to solve my problem? Thanks again!
I replaced "eType=gdal.GDT_UInt16" to "eType=gdal.GDT_Float32" and it worked. My binmask is an array of int32. Thank you.
Sorry for asking a question again. My clipped raster layer does not maintain the original color ramp of the input raster layer(ex. original raster: 0-2500, clipped raster: 4-40). Do you know how can I save a clipped raster to contain the original color ramp?
Just store the min/max values of your original raster and use these as limits for your clipped raster's color ramp.
Hello,
I have a question, How can we georeference a raster by corner points?
How helpful video ! thx a lot. Do you have any idea why my .tif plot doesn't show the color by the elevation as yours but just solid color?
Did you check the min/max values of your array? Could be that you have some very low nodata values (e.g. -9999). If that is the case, you need to adjust the colorbar in order to see something. Or just replace the low values with np.nan.
@@makingsenseremotely6207 I adjusted the colorbar and it's working now ! Thank you so much :)
Is there any way to get access to the files you use here? I want to code along but don't have a tif that works like this one. Thanks!
I have a .grc file, is there any way to use GDAL for this kind of raster file?
I'm trying to follow this tutorial on my Windows PC, but when I run this code on spyder, I get the "ModuleNotFoundError: No module named 'osgeo' " error. But I already have osgeo installed in my pc.
What do you mean by "you have osgeo installed"? You got the OSGeo4W Network Installer? You should be able to run Python and import gdal in the osgeo shell then. If you want to run this code in spyder: What worked for me on Windows was to download the GDAL wheel package. Check my GDAL Introduction + Installation video for more info.
How to get tif file
my python version is >3.8 and its returning error to import gdal from osgeo
what should I do
How to install osgeo ?
You need to install the GDAL Python bindings, which include the modules gdal, ogr, osr, gdalnumeric and gdalconst. Have a look at my gdal introduction + installation video.
Thank you soo much