“ GDAL processes raster data.”
Here, I will share some content on processing raster data with GDAL.
Sample data can still be downloaded from the following link:
https://pan.baidu.com/s/1bAnUo0S_ojxXdkyBqWAnLg?pwd=gxu2 Extraction code: gxu2
01
—
Data Reading
When importing the GDAL module, Python automatically registers all known GDAL drivers to support reading various compatible formats. Common file formats include TIFF, GeoTIFF, ASCII Grid, and Erdas Imagine’s img format.
We will open a band from the example Landsat8 data.
from osgeo import gdal# Adjust the path according to the local file storage locationfilepath = ".../Landsat8_sr_band4.tif"raster = gdal.Open(filepath)type(raster)osgeo.gdal.DatasetAt this point, the satellite image is stored as a GDAL Dataset object in the variable raster. We can check various attributes of the file.
raster.GetProjection()PROJCS["WGS 84 / UTM zone 35N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator",PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32635"]]raster.RasterXSize1045raster.RasterYSize1044raster.RasterCount1raster.GetMetadata(){'Band_1': 'band 4 surface reflectance', 'AREA_OR_POINT': 'Area'}We will store the GDAL raster band object in the variable band, and we can also check the properties of the band.
band = raster.GetRasterBand(1)type(band)osgeo.gdal.Bandgdal.GetDataTypeName(band.DataType)'Int16'We can calculate the statistics of the raster.
if band.GetMinimum() is None or band.GetMaximum()is None: band.ComputeStatistics(0) print("Statistics have been calculated.")band.GetMetadata(){'STATISTICS_MINIMUM': '-465','STATISTICS_MAXIMUM': '6346', 'STATISTICS_MEAN': '434.39829877723', 'STATISTICS_STDDEV': '381.09405826776','STATISTICS_VALID_PERCENT': '100'}print (band.GetNoDataValue())-9999.0print (band.GetMinimum())-465.0print (band.GetMaximum())6346GDAL is a powerful library for processing geospatial raster data, but it does not provide many functions for calculations. To perform more advanced calculations, we will read the raster data into a numerical array to leverage the capabilities of the NumPy library.
We can use the ReadAsArray() function to convert an existing GDAL Dataset or Band into a NumPy array.
rasterArray = raster.ReadAsArray()The GDAL library also comes with a gdal_array module, which acts as an interface between NumPy and GDAL. gdal_array can directly read and write raster files as NumPy arrays from files.
from osgeo import gdal_arrayrasterArray = gdal_array.LoadFile(filepath)02
—
GDAL Command Line Tools (Clipping, Stacking, and Band Calculations)
GDAL has powerful data conversion and processing capabilities, but these functions are not implemented directly in the library. For example, gdalwarp for clipping, mosaicking, and reprojection, gdal_merge.py for image mosaicking and stacking, and gdal_calc.py for raster calculations. These tools need to be run from the terminal/command prompt or as subprocesses in Python.
gdalwarp is a very useful tool for quickly clipping images, and we can use the -te option to specify the desired extent of the output file.
gdalwarp -te xmin ymin xmax ymax inputfile.tif outputfile.tifWe will repeat the clipping operation for all bands in the example data at once. We will generate corresponding commands for each band.
import globimport osFileList = glob.glob(os.path.join('.../Landsat8','*band*.tif'))# Replace with the local installation path of gdalwarp.gdal_path = ".../anaconda3/Library/bin/gdalwarp.exe"xmin, ymin, xmax, ymax = (472000, 6440000, 490000, 6455000)command = ""for fp in FileList: inputfile = fp outputfile = inputfile[:-4] + "_clip.tif" command += f'"{gdal_path}" -te {xmin} {ymin} {xmax} {ymax} "{inputfile}" "{outputfile}" \n'cmd_file = "ClipTurkufromLandsat.bat"f = open(os.path.join(cmd_file), 'w')f.write(command)f.close()Here’s a small issue: how to determine xmin, xmax, ymin, ymax. If we have the coordinates or vector of the desired clipping area, we can fill them directly or based on the vector coordinates. But now we are using an example, and we need to determine a value within the spatial extent of the raster data. We need to know the original raster’s xmin, xmax, ymin, ymax. In GDAL, we can calculate this based on the properties of the raster dataset (the raster variable mentioned above).
We first check the spatial properties using GetGeoTransform.
geo = raster.GetGeoTransformprint(geo)(468015.0, 30.0, 0.0, 6467415.0, 0.0, -30.0)The six numbers obtained represent the x-coordinate of the upper left corner (i.e., xmin), the horizontal resolution (pixel size in the x direction), the rotation parameter (usually 0), the y-coordinate of the upper left corner (i.e., ymax), the rotation parameter (usually 0), and the vertical resolution (pixel size in the y direction, usually negative, as geographic coordinates decrease downwards). Therefore, xmin, xmax, ymin, and ymax can be obtained using the following formulas:
xmin = geo[0]
ymax = geo[3]
xmax = geo[0] + width*geo[1]
ymin = geo[3] + height*geo[5]
where width is the value of RasterXSize, and height is the value of RasterYSize.
After running the above code, a file named ClipTurkufromLandsat.bat will be generated in the default working directory of our IDE or Jupyter notebook.
Then, we can run that file in the terminal window, or simply find the file and double-click to run it.
ClipTurkufromLandsat.batAfter clipping the image, we can stack bands 3 (green band), 4 (red band), and 5 (near-infrared band) to produce a false-color composite image. We can use the gdal_merge.py tool to merge these layers and specify whether we want to save each input layer as a separate band in the output file using the -separate option.
This time, we will try to run the command in Python as a subprocess.
import subprocess# Adjust the path according to the local location of python and gdal_merge, I use anaconda.python_path = r".../anaconda3/python.exe"gdal_merge_path = r".../anaconda3/Scripts/gdal_merge.py"# Adjust the path according to the local file storage locationinput_files = [ r".../landsat8/Landsat8_sr_band3_clip.tif", r".../landsat8/Landsat8_sr_band4_clip.tif", r".../landsat8/Landsat8_sr_band5_clip.tif"]output_file = r".../landsat8/Landsat8_GreenRedNir.tif"work_dir = r".../landsat8"command = [ python_path, gdal_merge_path, "-separate", "-o", output_file, *input_files]try: result = subprocess.run( command, cwd=work_dir, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True ) print("Merge successful! Output file:", output_file)except subprocess.CalledProcessError as e: print(f"Merge failed, error message:
{e.stderr}")except FileNotFoundError: print("Error: python or gdal_merge.py not found, please check the path.")<span>We can use gdal_calc.py</span> to perform simple repetitive calculations on raster data (e.g., band operations, statistical calculations, etc.). The basic syntax of gdal_calc.py is as follows.
gdal_calc.py -A input1.tif -B input2.tif [other_options] --outfile=outputfile.tifWe will run this command in Python as a subprocess to calculate NDVI.
# Adjust the path according to the local file storage locationred_band = r".../landsat8/Landsat8_sr_band4.tif"nir_band = r".../landsat8/Landsat8_sr_band5.tif"output_ndvi = r".../landsat8/Landsat8_NDVI.tif"gdal_calc_path = r".../anaconda3/Scripts/gdal_calc.py"command = [ "python", gdal_calc_path, "-A", nir_band, "-B", red_band, "--outfile", output_ndvi, "--calc", "(A - B) / (A + B)", "--NoDataValue", "-9999", "--type", "Float32"]try: subprocess.run(command, check=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) print(f"NDVI calculation completed, output file: {output_ndvi}")except subprocess.CalledProcessError as e: print(f"Calculation failed: {e.stderr}")If the raster files are large, you can use –co=”COMPRESS=LZW” parameter to reduce the output file size.