arcpy.NumPyArrayToRaster cont.

Question

I am still having difficulties converting a numpy array to a raster. 
I have created a raster for my cell sizes and have found the latitude min and longitude min for the lower_left_corner and use the following code:

from numpy import *
myarray=genfromtxt("E:/test/rasterise/myfile.txt")
mylon=genfromtxt("E:/test/rasterise/LATout.txt")
mylat=genfromtxt("E:/test/rasterise/LONout.txt")

corner=arcpy.Point(139.8, -39.2)
myRaster=arcpy.NumPyArrayToRaster(myarray, corner, mylon, mylat)

but I receive an error message:
Runtime error <type 'exceptions.TypeError'>: Argument x_cell_size: Expected a numeric value, Raster instance or path name.

however the documentation says that you can use a raster for the cell size. Not sure where I am going wrong here. 
Any feedback will be greatly appreciated.

SarahBurns · Answer

excellent thanks, that helps alot!

ChrisBater · Answer

This might help get you started with gdal....

Based on code I found here: http://www.gis.usu.edu/~chrisg/python/

'''
---------------------------------------------------------------------------
---------------------------------------------------------------------------
program: TOA to NDVI.py

Chris Bater 
GIS and Remote Sensing Technologist

Wildfire Resource Information Unit 
Wildfire Management Branch 
Forestry Division 
Alberta Sustainable Resource Development 
9th floor Great West Life Building 
9920 - 108 Street 
Edmonton, AB 
T5K 2M4 
Ph: 780-422-0669 
Cell: 780-446-1899
Email: chris.bater@gov.ab.ca

Date: September 2011

Purpose: Converts Landsat-5 TOA reflectance data to NDVI

---------------------------------------------------------------------------
----------------------------------------------------------------------------
'''

#----------------------------------------------------------------------
#Identify the directory containing the landsat data
workspace = r"C:\_LOCALdata	emp"
#----------------------------------------------------------------------

#----------------------------------------------------------------------
#Import and load modules
#----------------------------------------------------------------------

# Import system modules
import os, sys
from osgeo import gdal
import numpy as np
import numpy.ma as ma

#register the GDAL drivers
gdal.AllRegister()
driver = gdal.GetDriverByName('GTiff') #not sure if this is necessary
driver.Register()#not sure if this is necessary

print ""
print "All systems go"
print "Hit CTRL + c at any time to stop processing"
print ""

#----------------------------------------------------------------------
#Iterate through the landsat scenes and convert to TOA
#----------------------------------------------------------------------

try:
    
    files = os.listdir(workspace)
    for f in files:
        if not f.endswith('tif'):
            continue
        name = f[:-4]
        
        #open the landsat TOA file
        TOA = gdal.Open(workspace + "/" + f)

#get the image size
        rows = TOA.RasterYSize
        cols = TOA.RasterXSize
        bands = TOA.RasterCount
        
        if bands != 6:
            continue
        print "processing " + f
        
        #create an empty tif file
        outData = driver.Create(os.path.join(workspace, name + "_NDVI.tif"), cols, rows, 1, gdal.GDT_Float32)
        outBand = outData.GetRasterBand(1)

#assign the bands to variables
        b3 = TOA.GetRasterBand(3)
        b4 = TOA.GetRasterBand(4)
        
        
        #set the initial blocksize for analysis
        xBlockSize = 5000
        yBlockSize = 5000

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

#read bands into arrays
                d3 = b3.ReadAsArray(j, i, numCols, numRows).astype(np.float32)
                d4 = b4.ReadAsArray(j, i, numCols, numRows).astype(np.float32)
                
                #mask no data values (no data = 0)
                m3 = ma.masked_equal(d3, 0)
                m4 = ma.masked_equal(d4, 0)              
                del d3, d4
                
                ndvi = (m4 - m3) /(m4 + m3)
                del m3, m4
        
                outBand.WriteArray(ndvi, j, i)

# flush data to disk, set the NoData value and calculate stats
        outBand.FlushCache()
        outBand.SetNoDataValue(0)
        stats = outBand.GetStatistics(0, 1)

# georeference the image and set the projection
        outData.SetGeoTransform(TOA.GetGeoTransform())
        outData.SetProjection(TOA.GetProjection())

#memory management
        del TOA, outData, outBand

print 'script complete'
        
        
        
except Exception, e:

# If an error occurred, print line number and error message
    tb = sys.exc_info()[2]
    print "Line %i" % tb.tb_lineno
    print e.message

SarahBurns · Answer

It sounds like I will need to spend some time working out how to do this in GDAL....
thanks for the replies!