Garder la référence spatiale en utilisant arcpy.RasterToNumPyArray?

13

J'utilise ArcGIS 10.1 et je souhaite créer un nouveau raster basé sur deux rasters préexistants. Le RasterToNumPyArray a un bon exemple que je veux adapter.

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

Le problème est qu'il supprime la référence spatiale et également la taille des cellules. J'ai pensé que cela devait faire arcpy.env, mais comment les définir en fonction du raster en entrée? Je ne peux pas le comprendre.


En prenant la réponse de Luke, ceci est ma solution provisoire.

Les deux solutions de Luke définissent correctement la référence spatiale, l'étendue et la taille des cellules. Mais la première méthode ne portait pas correctement les données dans le tableau et le raster en sortie est rempli de nodata partout. Sa deuxième méthode fonctionne principalement, mais là où j'ai une grande région de nodata, elle se remplit de zéros en bloc et de 255. Cela peut avoir à voir avec la façon dont j'ai géré les cellules nodata, et je ne sais pas trop comment je le faisais (cela devrait être un autre Q cependant). J'ai inclus des images de ce dont je parle.

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

Image montrant les résultats. Dans les deux cas, les cellules nodales sont jaunes.

La deuxième méthode de Luke La deuxième méthode de Luke

Ma méthode provisoire Ma méthode provisoire

yosukesabai
la source

Réponses:

15

Découvrez la méthode Describe .

Quelque chose comme ce qui suit devrait fonctionner.

#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: La méthode arcpy.NumPyArrayToRaster prend un paramètre value_to_nodata. Utilisez-le comme ceci:

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
la source
Luke, Merci pour la réponse. Il me semble que je dois faire quelque chose d'hybride de vos deux méthodes? Les deux méthodes définissent l'étendue, la spref, etc., mais ne parviennent pas à mettre correctement les données en tableau ... La première méthode fait que toutes les cellules sont des nodata. La deuxième méthode fonctionne mieux, mais ils ne semblent pas gérer correctement les cellules nodata dans myArray (désolé, je n'ai pas dit que ma cellule contient des nodata et je veux les conserver intactes). Certains essais et erreurs m'ont fait penser que je dois adopter une deuxième approche, mais arcpy.env.outCoordinateSystem donne une apparence raisonnable à la sortie. pas beaucoup de comprendre comment.
yosukesabai
1
Vous n'avez pas posé de questions sur NoData, vous avez posé des questions sur la référence spatiale et la taille des cellules. La méthode arcpy.NumPyArrayToRaster prend un paramètre value_to_nodata.
user2856
Luke, Merci d'avoir édité. J'ai essayé votre méthode pour fournir le 5ème argument (value_to_nodata), mais j'ai quand même donné le chiffre en haut (cellule nodata remplie de 0 ou 255, et nodata_value n'est pas définie pour le raster en sortie). la seule solution de contournement que j'ai trouvée était de définir env.outputCoordinateSystem avant NumPyArrayToRaster, au lieu d'utiliser DefineProjection_management par la suite. Cela n'a pas de sens pourquoi cela fonctionne, mais je vais juste avec ma solution. Merci pour votre aide.
yosukesabai
1

J'ai rencontré des problèmes pour que ArcGIS gère correctement les valeurs NoData avec les exemples présentés ici. J'ai étendu l'exemple de blog reomtesensing.io (qui est plus ou moins similaire aux solutions présentées ici) pour mieux gérer NoData.

Apparemment, ArcGIS (10.1) aime la valeur -3.40282347e + 38 comme NoData. Je convertis donc en va-et-vient entre numpy NaN et -3.40282347e + 38. Le code est ici:

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
la source