Extracting Longitude into new raster

AbdelrazekElnashar — Mon, 29 Jun 2015 20:43:44 GMT

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

Re: Extracting Longitude into new raster

DarrenWiens2 — Mon, 29 Jun 2015 20:50:17 GMT

Since this is not ArcGIS-specific, you'll have better luck getting an answer on GIS StackExchange.

I will also say, when you post this over there, you should expand on your problem. I'm not sure what your question is.

Re: Extracting Longitude into new raster

XanderBakker — Mon, 29 Jun 2015 21:30:43 GMT

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

topic Extracting Longitude into new raster in Spatial Data Science Questions

Extracting Longitude into new raster

Re: Extracting Longitude into new raster

Re: Extracting Longitude into new raster