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