How to open, edit and save raster files (.tif) using Python
Fellow researchers and open-source GIS enthusiasts,
Welcome to my blog!
I’d like to start with a disclaimer – I may be a researcher of this very area but that doesn’t mean everything I do or write here will work for you, in your own desktop configurations and package versions. I have no responsibility if you lose data or mess up your installation. I also do not authorize any copies of my content.
Today, I am going to write about how to open, edit and save raster files using Python if you are familiarized with Python (mainly Anaconda).
Maybe you need to run every pixel in a raster through a Python function, or you want to do custom operations with them that are not available on GUI software such as QGIS. It also can be used to do batch processing. I showed another way to do batch processing here.
We are going to use GDAL on Python, installed through Conda, to do this.
How to do it
First thing you will need to do: install GDAL on Python using Anaconda
conda install gdal
I read it can be installed using pip as well, but have not tried it:
pip install gdal
and then load the package gdal from osgeo on your Python console or script. I use Anaconda Spyder.
from osgeo import gdal
Knowing the location of your raster to be opened, you can open it using the command:
step1 = gdal.Open('path_of_the_file.tif', gdal.GA_ReadOnly)
for opening the raster read-only and saving it on step1 variable.
Copy the transformation to a variable, it will be useful later.
GT_input = step1.GetGeoTransform()
Read the bands of your raster using GetRasterBand
step2 = step1.GetRasterBand(1)
Transform the read band into an array, to work with it.
img_as_array = step2.ReadAsArray()
Read the size of your array
Create the output to be saved with the same dimensions
import numpy as np output=np.zeros(shape=(size1,size2))
Do the processing you need to and save it to the output variable. Example, the whole raster to the power of 1.2:
for i in range(0,size1): for j in range(0,size2): output[i,j]=img_as_array[i,j] ** 1.2
Then, we will export the output variable into a file using rasterio.
You may need to install rasterio by using:
conda install -c conda-forge rasterio
Create a variable with the projection of your raster. Example:
Transform your output to a float
output = np.float32(output)
And finally, write the variable output to your file
with rasterio.open( 'output_map.tif', 'w', driver='GTiff', height=output.shape, width=output.shape, count=1, dtype=np.float32, crs=dst_crs, transform=GT_input, ) as dest_file: dest_file.write(output, 1) dest_file.close()
This example saves an output called output_map.tif
Example with a random raster
I generated a random raster by binomial distribution on QGIS 3.18, with values between 0 and 10. The actual range of the generated raster was [2,7].
After running the code on Python, the following raster was generated:
It looks the same, but the range has changed:
Showing that the code worked!
Edit (August 2nd, 2021): if you are having issues with GDAL-style transforms in rasterio, there is a quick solve for it. Check out this post.
- How to read grib2 raster data on Python
- Solved: new version rasterio 1.0+ won’t let me use GDAL transforms
- Batch processing – automatizing the use of GDAL and SAGA GIS tools using Bash Shell scripts
- How to load (many) raster layers on QGIS 3.18 by batch process and alter their symbology together
- Solved: QGIS 3.20 Python Console shows an error when running Clip Vector by Extent tool