Raumbezug mit arcpy.RasterToNumPyArray beibehalten?

13

Ich verwende ArcGIS 10.1 und möchte ein neues Raster erstellen, das auf zwei bereits vorhandenen Rastern basiert. Das RasterToNumPyArray hat ein gutes Beispiel, das ich anpassen möchte.

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 ist, dass es den Raumbezug und auch die Zellengröße entfernt. Ich dachte, es muss arcpy.env sein, aber wie stelle ich sie basierend auf dem Eingabe-Raster ein? Ich kann es nicht herausfinden.


Nach Lukes Antwort ist dies meine vorläufige Lösung.

Beide von Lukes Lösung eingestellten Raumbezug, Ausdehnung und Zellgröße korrekt. Bei der ersten Methode wurden die Daten im Array nicht korrekt übertragen, und das Ausgabe-Raster ist überall mit Knotendaten gefüllt. Seine zweite Methode funktioniert meistens, aber wo ich große Region der Knoten habe, füllt es mit blocky Nullen und 255s. Dies kann damit zusammenhängen, wie ich mit Knotenzellen umgegangen bin, und ich bin mir nicht ganz sicher, wie ich es gemacht habe (sollte aber ein anderes Q sein). Ich habe Bilder von dem beigefügt, wovon ich spreche.

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

Bild zeigt Ergebnisse. In beiden Fällen sind die Knotenzellen gelb dargestellt.

Lukes zweite Methode Lukes zweite Methode

Meine vorläufige Methode Meine vorläufige Methode

Yosukesabai
quelle

Antworten:

15

Schauen Sie sich die Describe- Methode an.

Etwas wie das Folgende sollte funktionieren.

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

ODER

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

EDIT: Die Methode arcpy.NumPyArrayToRaster akzeptiert einen value_to_nodata-Parameter. Benutze es so:

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)
user2856
quelle
Luke, danke für die Antwort. Es scheint mir, dass ich etwas Hybrides aus Ihren beiden Methoden machen muss? Beide Methoden setzen Extent, SPREF usw. korrekt, aber ordnen die Daten nicht korrekt an ... Die erste Methode macht alle Zellen zu Knotendaten. Die zweite Methode funktioniert besser, aber sie scheinen nicht mit Knotenzellen in myArray richtig umzugehen (sorry, dass ich nicht gesagt habe, dass meine Zelle Knotendaten hat und ich diese intakt halten möchte). Einige Versuche ließen mich denken, ich müsse einen zweiten Ansatz wählen, aber arcpy.env.outCoordinateSystem lässt die Ausgabe vernünftig erscheinen. obwohl nicht viel Verständnis dafür, wie so.
Yosukesabai
1
Sie haben nicht nach NoData gefragt, sondern nach Raumbezug und Zellengröße. Die arcpy.NumPyArrayToRaster- Methode akzeptiert einen value_to_nodata-Parameter.
User2856
Luke, vielen Dank für die Bearbeitung. Ich habe versucht, mit Ihrer Methode das fünfte Argument (value_to_nodata) bereitzustellen, aber ich habe trotzdem die Zahl oben erhalten (Knotenzelle gefüllt mit 0 oder 255, und nodata_value ist nicht für das Ausgabe-Raster festgelegt). Die einzige Problemumgehung, die ich gefunden habe, war das Festlegen von env.outputCoordinateSystem vor NumPyArrayToRaster, anstatt danach DefineProjection_management zu verwenden. Es macht keinen Sinn, warum dies funktioniert, aber ich gehe einfach mit meiner Lösung. Danke für all die Hilfe.
Yosukesabai
1

Ich hatte einige Probleme, ArcGIS dazu zu bringen, NoData-Werte mit den hier gezeigten Beispielen korrekt zu verarbeiten. Ich habe das Beispiel aus dem Blog reomtesensing.io (das den hier gezeigten Lösungen mehr oder weniger ähnlich ist) erweitert, um besser mit NoData umgehen zu können.

Anscheinend gefällt ArcGIS (10.1) der Wert -3.40282347e + 38 als NoData. Also konvertiere ich zwischen numpy NaN und -3.40282347e + 38 hin und her. Der Code ist hier:

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)
Mathiaskopo
quelle