I am using following code to create bt2.vrt file using python.
import numpy as np import gdal from gdalconst import * import osr # Function to read the original file's projection: def GetGeoInfo(FileName): SourceDS = gdal.Open(FileName, GA_ReadOnly) xsize = SourceDS.RasterXSize ysize = SourceDS.RasterYSize GeoT = SourceDS.GetGeoTransform() Projection = osr.SpatialReference() Projection.ImportFromWkt(SourceDS.GetProjectionRef()) DataType = SourceDS.GetRasterBand(1).DataType #DataType = gdal.GetDataTypeName(DataType) return xsize, ysize, GeoT, Projection, DataType # Function to write a new file. def Create_vrt(Name, Array, driver,xsize, ysize, GeoT, Projection, DataType): if DataType == 'UInt16': DataType = gdal.GDT_Float32 NewFileName = 'bt2'+'.vrt' # Set nans to the original No Data Value # Set up the dataset DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType ) # the '1' is for band 1. DataSet.SetGeoTransform(GeoT) DataSet.SetProjection( Projection.ExportToWkt() ) # Write the array DataSet.GetRasterBand(1).WriteArray( Array ) return NewFileName # Open the original file FileName='/home/vigna/TIFFs/bt.tif' DataSet = gdal.Open(FileName, GA_ReadOnly) # Get the first (and only) band. Band = DataSet.GetRasterBand(1) # Open as an array. Array = Band.ReadAsArray() # Now I do some processing on Array, it's pretty complex NewArray = Array print(NewArray) #save the new file, in the meantime I have xsize, ysize, GeoT, Projection, DataType = GetGeoInfo(FileName) # Set up the vrt driver driver = gdal.GetDriverByName('VRT') # Now turn the array into a vrt. NewFileName = Create_vrt('bt3', NewArray, driver, xsize, ysize, GeoT, Projection, DataType) #Create_vrt() print("ok")
bt2.vrt file created using above code is : : https://i.stack.imgur.com/JYqVf.png
According to : Unable to warp HDF5 files I need to include lat.vrt and long.vrt as y_dataset and X_dataset respectively to perform "gdal_warp" on bt.vrt file. How to include it using python code?