Extracting Longitude into new raster

AbdelrazekElnashar

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")

XanderBakker

If you are willing to do this within ArcGIS (and if I understand your question correctly) you could have a look at this post NumPyArray to Raster And Back Again by Curtis Price