Utrzymywanie odniesienia przestrzennego za pomocą arcpy.RasterToNumPyArray?

Korzystam z ArcGIS 10.1 i chcę utworzyć nowy raster na podstawie dwóch istniejących rastrów. RasterToNumPyArray ma dobry przykład, który chcę, aby dostosować.

import arcpy
import numpy 
myArray = arcpy.RasterToNumPyArray('C:/data/inRaster')
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum
newRaster = arcpy.NumPyArrayToRaster(myArrayPerc)
newRaster.save("C:/output/fgdb.gdb/PercentRaster")

Problem polega na tym, że usuwa odniesienie przestrzenne, a także rozmiar komórki. Pomyślałem, że to musi zrobić arcpy.env, ale jak ustawić je na podstawie rastra wejściowego? Nie mogę tego rozgryźć.

Biorąc odpowiedź Luke'a, to jest moje wstępne rozwiązanie.

Oba rozwiązania Luke'a poprawnie ustawiły odniesienie przestrzenne, zasięg i rozmiar komórki. Ale pierwsza metoda nie przeprowadziła poprawnie danych w tablicy, a wyjściowy raster jest wszędzie zapełniony nodanymi. Jego druga metoda działa głównie, ale tam, gdzie mam duży obszar nodata, wypełnia się blokowymi zerami i 255s. Może to mieć związek z tym, jak obsłużyłem komórki nodata, i nie jestem do końca pewien, jak to robiłem (powinno to być kolejne Q). Dołączyłem zdjęcia tego, o czym mówię.

#Setting the raster properties directly 
import arcpy 
import numpy 

inRaster0='C:/workspace/test0.tif' 
inRaster1='C:/workspace/test1.tif' 
outRaster='C:/workspace/test2.tif' 

dsc=arcpy.Describe(inRaster0) 
sr=dsc.SpatialReference 
ext=dsc.Extent 
ll=arcpy.Point(ext.XMin,ext.YMin) 

# sorry that i modify calculation from my original Q.  
# This is what I really wanted to do, taking two uint8 rasters, calculate 
# the ratio, express the results as percentage and then save it as uint8 raster.
tmp = [ np.ma.masked_greater(arcpy.RasterToNumPyArray(_), 100) for _ in inRaster0, inRaster1]
tmp = [ np.ma.masked_array(_, dtype=np.float32) for _ in tmp]
tmp = ((tmp[1] ) / tmp[0] ) * 100
tmp = np.ma.array(tmp, dtype=np.uint8)
# i actually am not sure how to properly carry the nodata back to raster...  
# but that's another Q
tmp = np.ma.filled(tmp, 255)

# without this, nodata cell may be filled with zero or 255?
arcpy.env.outCoordinateSystem = sr

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight) 

newRaster.save(outRaster) 

Obraz przedstawiający wyniki. W obu przypadkach komórki nodata mają kolor żółty.

Druga metoda Luke'a

Moja wstępna metoda

Sprawdź metodę opisywania .

Powinno działać coś takiego.

#Using arcpy.env
import arcpy
import numpy

inRaster='C:/workspace/test1.tif'
outRaster='C:/workspace/test2.tif'

dsc=arcpy.Describe(inRaster)
arcpy.env.extent=dsc.Extent
arcpy.env.outputCoordinateSystem=dsc.SpatialReference
arcpy.env.cellSize=dsc.meanCellWidth

myArray = arcpy.RasterToNumPyArray(r)
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc)
newRaster.save(outRaster)

LUB

#Setting the raster properties directly
import arcpy
import numpy

inRaster='C:/workspace/test1.tif'
outRaster='C:/workspace/test2.tif'

dsc=arcpy.Describe(inRaster)
sr=dsc.SpatialReference
ext=dsc.Extent
ll=arcpy.Point(ext.XMin,ext.YMin)

myArray = arcpy.RasterToNumPyArray(inRaster)
myArraySum = myArray.sum(1)
myArraySum.shape = (myArray.shape[0],1)
myArrayPerc = (myArray * 1.0)/ myArraySum

newRaster = arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight)
arcpy.DefineProjection_management(newRaster, sr)

newRaster.save(outRaster)

EDYCJA: Metoda arcpy.NumPyArrayToRaster przyjmuje parametr value_to_nodata. Używaj go w ten sposób:

try:
    noDataValue=dsc.noDataValue
    arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight,noDataValue)
except AttributeError: #no data is not defined
    arcpy.NumPyArrayToRaster(myArrayPerc,ll,dsc.meanCellWidth,dsc.meanCellHeight)

Miałem pewne problemy z tym, że ArcGIS poprawnie obsługiwał wartości NoData z przykładami pokazanymi tutaj. Rozszerzyłem przykład z bloga reomtesensing.io (który jest mniej więcej podobny do pokazanych tutaj rozwiązań), aby lepiej obsługiwać NoData.

Najwyraźniej ArcGIS (10.1) lubi wartość -3,40282347e + 38 jako NoData. Więc przeliczam tam iz powrotem między numpy NaN a -3,40282347e + 38. Kod jest tutaj:

import arcpy
import numpy as np

infile  = r'C:\data.gdb\test_in'
outfile = r'C:\data.gdb\test_out'

# read raster with No Data as numpy NaN
in_arr  = arcpy.RasterToNumPyArray(infile,nodata_to_value = np.nan)

# processing
out_arr = in_arr * 10

# convert numpy NaN to -3.40282347e+38
out_arr[np.isnan(out_arr)] = -3.40282347e+38

# information on input raster
spatialref = arcpy.Describe(infile).spatialReference
cellsize1  = arcpy.Describe(infile).meanCellHeight
cellsize2  = arcpy.Describe(infile).meanCellWidth
extent     = arcpy.Describe(infile).Extent
pnt        = arcpy.Point(extent.XMin,extent.YMin)

# save raster
out_ras = arcpy.NumPyArrayToRaster(out_arr,pnt,cellsize1,cellsize2, -3.40282347e+38)
out_ras.save(outfile)
arcpy.DefineProjection_management(outfile, spatialref)