GDAL Usage

---

GDAL write GEOtiff

from osgeo import gdal
from osgeo import osr
import os
import numpy as np

def write_geo_tiff(img, path, scale, dataType = 'Float32',NoData=-999,BandNames=None):
    
    img[np.isnan(img)] = NoData
    if dataType == 'Float32':
        dType = gdal.GDT_Float32
        img = img.astype(np.float32)
        
    elif dataType == 'Int32':
        dType = gdal.GDT_Int32
        img = img.astype(np.int32)
    
    elif dataType == 'Int16':
        dType = gdal.GDT_Int16
        img = img.astype(np.int16)
    
    elif dataType == 'Uint8':
        dType = gdal.GDT_Byte
        img = img.astype(np.uint8)
        
    else:
        raise ValueError('type error')
        
    # create transformation parameters
    lon = img[-2,:,:]
    lat = img[-1,:,:]
    upLeft_x = np.min(lon)
    scale_h = scale
    rotation = 0
    upLeft_y = np.max(lat)
    scale_v = -1*scale
    im_geotrans = [upLeft_x,scale_h,rotation,upLeft_y,rotation,scale_v]  
    
    # create crs
    crs = osr.SpatialReference()
    crs.ImportFromEPSG(4326)
    
    # creat geotif
    if len(img.shape) == 2:
        img = img[np.newaxis,:,:]
    driver = gdal.GetDriverByName("GTiff")
 
    dataset = driver.Create(path, 
                            img.shape[2], img.shape[1], img.shape[0], dType)
    if(dataset!= None):
        dataset.SetGeoTransform(im_geotrans)
        dataset.SetProjection(crs.ExportToWkt())
    for i in range(img.shape[0]):
        band = dataset.GetRasterBand(i+1)
        band.SetNoDataValue(NoData)
        if not BandNames is None:
            band.SetDescription('Band%d_%s'%(i+1,BandNames[i]))
        band.WriteArray(img[i,:,:])
    del dataset

Clip

def gdal_clip(input_shape, inRaster, out_raster):

	inRaster=gdal.Open(inRaster)
	# execute  
	ds = gdal.Warp(out_raster,
		  inRaster,
		  format = 'GTiff',
		  cutlineDSName = input_shape,      # or any other file format
		  # cutlineWhere="FIELD = 'whatever'",# commented out by Qi, optionally you can filter your cutline (shapefile) based on attribute values
		  cropToCutline = True, # added by Qi, 2020/07/26
		  dstNodata = 0)              # select the no data value you like
	ds = None     #do other stuff with ds object, it is your cropped dataset. in this case we only close the dataset.

Back