Mantendo referência espacial usando arcpy.RasterToNumPyArray?

13

Estou usando o ArcGIS 10.1 e quero criar uma nova varredura com base em duas rasters existentes. O RasterToNumPyArray tem um bom exemplo que eu quero adaptar.

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

O problema é que retira a referência espacial e também o tamanho da célula. Achei que tinha que fazer arcpy.env, mas como faço para configurá-los com base na entrada raster? Eu não consigo entender.

Tomando a resposta de Luke, esta é minha solução provisória.

As duas soluções de Lucas definiram a referência espacial, a extensão e o tamanho da célula corretamente. Mas o primeiro método não carregava dados na matriz corretamente e a varredura de saída é preenchida com nodata em todos os lugares. Seu segundo método funciona principalmente, mas onde eu tenho uma grande região de nodata, ele preenche com zeros em bloco e 255s. Isso pode ter a ver com a maneira como lidei com as células nodata, e não tenho muita certeza de como estava fazendo isso (deveria ser outro Q). Incluí imagens do que estou falando.

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

Imagem mostrando resultados. Nos dois casos, as células nodatas são mostradas em amarelo.

Segundo método de Lucas Segundo método de Lucas

Meu método experimental Meu método experimental

arcpy raster arcgis-10.1 coordinate-system numpy yosukesabai
fonte

Respostas:

15

Confira o método Descrever .

Algo como o seguinte deve funcionar.

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

OU

#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: O método arcpy.NumPyArrayToRaster usa um parâmetro value_to_nodata. Use-o assim:

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
fonte
Luke, obrigado pela resposta. Parece-me que preciso fazer algo híbrido dos seus dois métodos? Ambos os métodos definem extensão, spref etc. corretos, mas falham em organizar os dados corretamente ... O primeiro método faz com que toda a célula seja nodata. O segundo método funciona melhor, mas eles parecem não manipular corretamente as células nodata no myArray (desculpe-me por não ter dito que meu celular tem alguns nodata e quero mantê-los intactos). Algumas tentativas e erros me fizeram pensar que eu tinha que adotar a segunda abordagem, mas o arcpy.env.outCoordinateSystem faz com que a saída pareça razoável. não há muita compreensão de como isso acontece.
yosukesabai
1
Você não perguntou sobre NoData, você perguntou sobre referência espacial e tamanho da célula. O método arcpy.NumPyArrayToRaster usa um parâmetro value_to_nodata.
user2856
Luke, obrigado pela edição. Eu tentei o seu método de fornecer o 5º argumento (value_to_nodata), mas ainda assim produzi a figura na parte superior (célula nodata preenchida com 0 ou 255 e nodata_value não está definido para a varredura de saída). a única solução alternativa que encontrei foi definir env.outputCoordinateSystem antes de NumPyArrayToRaster, em vez de usar DefineProjection_management posteriormente. Não faz sentido porque isso funciona, mas eu apenas uso a minha solução. Obrigado por toda a ajuda.
yosukesabai
1

Eu tive alguns problemas para conseguir que o ArcGIS manipulasse os valores NoData corretamente com os exemplos mostrados aqui. Estendi o exemplo do blog reomtesensing.io (que é mais ou menos semelhante às soluções mostradas aqui) para lidar melhor com o NoData.

Aparentemente, o ArcGIS (10.1) gosta do valor -3.40282347e + 38 como NoData. Portanto, eu converto entre NaN e -3,40282347e + 38. O código está aqui:

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
fonte