topic Re: arcpy.NumPyArrayToRaster cont. in Python Questions

arcpy.NumPyArrayToRaster cont.

SarahBurns — Wed, 25 Jan 2012 02:53:43 GMT

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.

Re: arcpy.NumPyArrayToRaster cont.

curtvprice — Wed, 25 Jan 2012 04:50:41 GMT

however the documentation says that you can use a raster for the cell size. Not sure where I am going wrong here.

By "raster" I believe they mean an ArcGIS-recognized raster dataset or raster band, not a numpy array.

Re: arcpy.NumPyArrayToRaster cont.

ChrisBater — Wed, 25 Jan 2012 13:57:23 GMT

I've been struggling with this as well. My work-around has been to use gdal. However, if I'm forced to use arcpy, then I crank out the array as a raster with no spatial reference, then use "rescale," "shift," and "define projection" to put the final surface back to where it should be. I wrote a function to do this, but it's definitely not optimal. There needs to be a better, more in depth usage example included in the numpyarraytoraster tool's help doc.

Re: arcpy.NumPyArrayToRaster cont.

SarahBurns — Thu, 26 Jan 2012 23:31:39 GMT

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

Re: arcpy.NumPyArrayToRaster cont.

ChrisBater — Sat, 11 Dec 2021 07:09:48 GMT

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\temp"
#----------------------------------------------------------------------




#----------------------------------------------------------------------
#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

Re: arcpy.NumPyArrayToRaster cont.

SarahBurns — Tue, 31 Jan 2012 02:24:35 GMT

excellent thanks, that helps alot!