Menyimpan referensi spasial menggunakan arcpy.RasterToNumPyArray?

13

Saya menggunakan ArcGIS 10.1, dan ingin membuat raster baru berdasarkan dua raster yang sudah ada sebelumnya. The RasterToNumPyArray memiliki contoh yang baik yang saya ingin beradaptasi.

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

Masalahnya adalah ia menghilangkan referensi spasial dan juga ukuran sel. Saya pikir itu harus dilakukan arcpy.env, tetapi bagaimana cara mengaturnya berdasarkan input raster? Saya tidak bisa mengetahuinya.


Mengambil jawaban Luke, ini adalah solusi sementara saya.

Kedua solusi Luke mengatur referensi spasial, luas dan ukuran sel dengan benar. Tetapi metode pertama tidak membawa data dalam array dengan benar dan output raster diisi dengan nodata di mana-mana. Metode keduanya bekerja sebagian besar, tetapi di mana saya memiliki wilayah besar nodata, itu diisi dengan nol dan 255 gumpal. Ini mungkin ada hubungannya dengan bagaimana saya menangani sel-sel nodata, dan saya tidak yakin bagaimana saya melakukannya (seharusnya menjadi Q lain). Saya menyertakan gambar dari apa yang saya bicarakan.

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

Gambar menampilkan hasil. Saya berdua sel-sel nodata ditampilkan kuning.

Metode kedua Luke Metode kedua Luke

Metode tentatif saya Metode tentatif saya

yosukesabai
sumber

Jawaban:

15

Periksa metode Jelaskan .

Sesuatu seperti yang berikut ini seharusnya bekerja.

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

ATAU

#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: Metode arcpy.NumPyArrayToRaster mengambil parameter value_to_nodata. Gunakan seperti ini:

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)
pengguna2856
sumber
Luke, Terima kasih atas jawabannya. Tampak bagi saya bahwa saya perlu melakukan sesuatu hibrida dari dua metode Anda? Kedua metode menetapkan luas, spref dll. Benar, tetapi gagal untuk menyusun data dengan benar ... Metode pertama membuat semua sel menjadi nodata. Metode kedua bekerja lebih baik, tetapi mereka tampaknya tidak menangani sel-sel nodata di myArray dengan benar (maaf saya tidak memberi tahu sel saya memiliki beberapa nodata dan saya ingin menjaga sel-sel itu tetap utuh). Beberapa trial and error membuat saya berpikir bahwa saya harus mengambil pendekatan kedua tetapi arcpy.env.outCoordinateSystem membuat output terlihat masuk akal. tidak banyak pemahaman tentang bagaimana hal itu terjadi.
yosukesabai
1
Anda tidak bertanya tentang NoData, Anda bertanya tentang referensi spasial dan ukuran sel. The arcpy.NumPyArrayToRaster Metode mengambil parameter value_to_nodata.
user2856
Luke, Terima kasih sudah mengedit. Saya mencoba metode Anda memberikan argumen ke-5 (value_to_nodata), tetapi saya masih menghasilkan angka di atas (sel nodata diisi dengan 0 atau 255, dan nodata_value tidak diatur untuk raster keluaran). satu-satunya solusi yang saya temukan adalah mengatur env.outputCoordinateSystem sebelum NumPyArrayToRaster, alih-alih menggunakan DefineProjection_management sesudahnya. Tidak masuk akal mengapa ini bekerja, tetapi saya hanya pergi dengan solusi saya. terima kasih atas semua bantuannya.
yosukesabai
1

Saya punya beberapa masalah dalam mendapatkan ArcGIS untuk menangani nilai-nilai NoData dengan benar dengan contoh-contoh yang ditunjukkan di sini. Saya memperluas contoh dari blog reomtesensing.io (yang kurang lebih mirip dengan solusi yang ditampilkan di sini) untuk menangani NoData dengan lebih baik.

Rupanya ArcGIS (10.1) menyukai nilai -3.40282347e + 38 sebagai NoData. Jadi saya mengkonversi bolak-balik antara NaN numpy dan -3.40282347e + 38. Kode di sini:

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
sumber