Utrzymywanie odniesienia przestrzennego za pomocą arcpy.RasterToNumPyArray?

13

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 Druga metoda Luke'a

Moja wstępna metoda Moja wstępna metoda

yosukesabai
źródło

Odpowiedzi:

15

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)
użytkownik2856
źródło
Luke, dziękuję za odpowiedź. Wydaje mi się, że muszę zrobić coś hybrydowego z twoich dwóch metod? Obie metody ustawiają zasięg, spref itp. Są poprawne, ale nie układają poprawnie danych ... Pierwsza metoda powoduje, że wszystkie komórki są nodata. Druga metoda działa lepiej, ale wydaje się, że nie obsługują poprawnie komórek nodata w myArray (przepraszam, że nie powiedziałem, że moja komórka ma jakieś nodata i chcę je zachować). Niektóre próby i błędy sprawiły, że pomyślałem, że muszę zastosować drugie podejście, ale arcpy.env.outCoordinateSystem sprawia, że ​​dane wyjściowe wyglądają na rozsądne. choć nie bardzo rozumiem, jak to się dzieje.
yosukesabai
1
Nie pytałeś o NoData, pytałeś o odniesienie przestrzenne i rozmiar komórki. Metoda arcpy.NumPyArrayToRaster przyjmuje parametr value_to_nodata.
user2856
Luke, dziękuję za edycję. Wypróbowałem twoją metodę dostarczania piątego argumentu (wartość_do_nodaty), ale nadal otrzymałem cyfrę u góry (komórka nodata wypełniona 0 lub 255, a wartość nodata_value nie jest ustawiona dla wyjściowego rastra). jedynym obejściem, jakie znalazłem, było ustawienie env.outputCoordinateSystem przed NumPyArrayToRaster, zamiast używania DefineProjection_management później. Nie ma sensu, dlaczego to działa, ale po prostu idę z moim rozwiązaniem. dziękuję za całą pomoc.
yosukesabai
1

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)
Mathiaskopo
źródło