But the boundary not matched, what shall i do?
# -*- coding: utf-8 -*-
#!/usr/bin/env python
print("Started\n")
import gdal,osr
import numpy as np
#*******************************************************************************
inRaster=r"E:\SEBAL_METRIC_Approaches\ASD\gcsDEM.tif"
PrjRaster="E:\SEBAL_METRIC_Approaches\ASD\pcsDEM.tif"
outRaster=r"E:\SEBAL_METRIC_Approaches\Outputs1\Longfh.tif"
Driver="GTIFF" # HFA Erdas Imagine (.img), GTIFF (TIFF)
#*******************************************************************************
inRaster=gdal.Open(inRaster)
rows=inRaster.RasterYSize # Z.shape[0]
cols=inRaster.RasterXSize # Z.shape[1]
Z=inRaster.ReadAsArray()
gt=inRaster.GetGeoTransform()
inBand=inRaster.GetRasterBand(1)
(xBlockSize,yBlockSize)=inBand.GetBlockSize()
(xOrigin,pixelWidth,rotation,yOrigin,rotation,pixelHeight)=inRaster.GetGeoTransform()
X=np.array(np.arange(xOrigin,(xOrigin+(cols-0.5)*pixelWidth),pixelWidth))
Y=np.arange(yOrigin,(yOrigin+rows*pixelHeight),pixelHeight)
driver=gdal.GetDriverByName(Driver)
outRaster=driver.Create(outRaster,cols,rows,1,gdal.GDT_Float32)
outBand=outRaster.GetRasterBand(1)
for i in range(0,rows,yBlockSize):
if i+yBlockSize<rows:numRows=yBlockSize
else:numRows=rows-i
for j in range(0,cols,xBlockSize):
if j+xBlockSize<cols:numCols=xBlockSize
else:numCols=cols-j
BandArray=inBand.ReadAsArray(j,i,numCols,numRows).astype(np.float)
OnesArr=np.ones((rows,cols))
outBand.WriteArray(X,j,i)
outBand.FlushCache()
stats=outBand.GetStatistics(0,1)
PrjRaster=gdal.Open(PrjRaster)
(xOrigin,pixelWidth,rotation,yOrigin,rotation,pixelHeight)=PrjRaster.GetGeoTransform()
outRaster.SetGeoTransform((xOrigin,pixelWidth,rotation,yOrigin,rotation,pixelHeight))
outRaster.SetProjection(PrjRaster.GetProjection())
##gdal.SetConfigOption("HFA_USE_RRD","YES")
gdal.SetConfigOption("COMPRESS_OVERVIEW","LZW")
outRaster.BuildOverviews(overviewlist=[2,4,8,16,32,64,128])
print("\nFinished")
