Jak uzyskać współrzędne XY i wartość komórki każdego piksela w rastrze za pomocą Pythona?

Jestem naprawdę nowy w Pythonie i chciałbym wiedzieć, czy istnieje szybka metoda uzyskiwania wartości komórek piksela rastrowego piksel po pikselu i współrzędnych (współrzędna mapy XY środka każdego piksela) za pomocą Pythona w ArcGIS 10?

Aby opisać to dokładniej, muszę uzyskać mapę X, mapę Y i wartość komórki pierwszego piksela, przypisać te trzy wartości do trzech zmiennych i powtórzyć ten krok dla pozostałych pikseli (pętla przez cały raster).

Myślę, że muszę bardziej opisać moje pytanie. Problem w tym, że muszę uzyskać lokalizację XY piksela pierwszego rastra i uzyskać wartości komórek kilku innych rastrów odpowiadających tej lokalizacji XY. Proces ten powinien przebiegać w pętli przez każdy piksel pierwszego rastra bez tworzenia pliku kształtu pośredniego punktu, ponieważ zajmie to naprawdę dużo czasu, ponieważ muszę obsługiwać raster z prawie 8 miliardami pikseli. Ponadto muszę to zrobić za pomocą Pythona w ArcGIS 10.

@JamesS: Dziękuję bardzo za sugestię. Tak, to zadziałałoby dla jednego rastra, ale muszę również zebrać wartości komórek dla kilku innych rastrów. Problem polega na tym, że po uzyskaniu współrzędnych X i Y pierwszego piksela pierwszego rastra muszę uzyskać wartość komórki drugiego rastra odpowiadającą lokalizacji X, Y pierwszego rastra, a następnie trzeciego rastra i tak dalej. Tak więc myślę, że podczas zapętlania pierwszego rastra pobieranie lokalizacji X i Y piksela i uzyskiwanie wartości komórek drugiego rastra odpowiadających tej lokalizacji powinno odbywać się jednocześnie, ale nie jestem pewien. Można to zrobić, konwertując pierwszy raster na plik kształtu punktu i wykonując Wyodrębnij wielowartościowe do funkcji punktowej w ArcGIS 10, ale ja '

@hmfly: Dzięki, tak, ta metoda (RastertoNumpyarray) zadziała, jeśli uda mi się uzyskać współrzędne znanej wartości wiersza i kolumny tablicy.

@ whuber: Nie chcę wykonywać żadnych obliczeń, wszystko, co muszę zrobić, to zapisać współrzędne XY i wartości komórek w pliku tekstowym i to wszystko

Zgodnie z pomysłem @ Dango stworzyłem i przetestowałem (na małych rastrach o tym samym zasięgu i rozmiarze komórki) następujący kod:

import arcpy, numpy

inRaster = r"C:\tmp\RastersArray.gdb\InRaster"
inRaster2 = r"C:\tmp\RastersArray.gdb\InRaster2"

##Get properties of the input raster
inRasterDesc = arcpy.Describe(inRaster)

#coordinates of the lower left corner
rasXmin = inRasterDesc.Extent.Xmin
rasYmin = inRasterDesc.Extent.Ymin

# Cell size, raster size
rasMeanCellHeight = inRasterDesc.MeanCellHeight
rasMeanCellWidth = inRasterDesc.MeanCellWidth
rasHeight = inRasterDesc.Height
rasWidth = inRasterDesc.Width

##Calculate coordinates basing on raster properties
#create numpy array of coordinates of cell centroids
def rasCentrX(rasHeight, rasWidth):
    coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)
    return coordX
inRasterCoordX = numpy.fromfunction(rasCentrX, (rasHeight,rasWidth)) #numpy array of X coord

def rasCentrY(rasHeight, rasWidth):
    coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)
    return coordY
inRasterCoordY = numpy.fromfunction(rasCentrY, (rasHeight,rasWidth)) #numpy array of Y coord

#combine arrays of coordinates (although array for Y is before X, dstack produces [X, Y] pairs)
inRasterCoordinates = numpy.dstack((inRasterCoordY,inRasterCoordX))


##Raster conversion to NumPy Array
#create NumPy array from input rasters 
inRasterArrayTopLeft = arcpy.RasterToNumPyArray(inRaster)
inRasterArrayTopLeft2 = arcpy.RasterToNumPyArray(inRaster2)

#flip array upside down - then lower left corner cells has the same index as cells in coordinates array
inRasterArray = numpy.flipud(inRasterArrayTopLeft)
inRasterArray2 = numpy.flipud(inRasterArrayTopLeft2)


# combine coordinates and value
inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

#add values from second raster
rasterValuesArray = numpy.dstack((inRasterFullArray, inRasterArray2.T))

Na podstawie kodu @hmfly możesz uzyskać dostęp do pożądanych wartości:

(height, width, dim )=rasterValuesArray.shape
for row in range(0,height):
    for col in range(0,width):
        #now you have access to single array of values for one cell location

Niestety istnieje jedno „ale” - kod jest odpowiedni dla tablic NumPy, które mogą być obsługiwane przez pamięć systemową. W moim systemie (8 GB) największa tablica wynosiła około 9000,9000.

Ponieważ moje doświadczenie nie pozwala mi zapewnić więcej pomocy, możesz rozważyć kilka sugestii dotyczących postępowania z dużymi tablicami: /programming/1053928/python-numpy-very-large-matrices

arcpy.RasterToNumPyArrayMetoda pozwala określić podzbiór rastra przekonwertowanego na tablicę NumPy ( strona pomocy ArcGIS10 ), co może być przydatne, gdy fragmentuje się duży zestaw danych w podmacierze.

Jeśli chcesz tylko uzyskać wartości w pikselach (wiersz, kolumna), możesz napisać arkadowy skrypt w następujący sposób:

import arcpy
raster = arcpy.Raster("yourfilepath")
array = arcpy.RasterToNumPyArray(raster)
(height, width)=array.shape
for row in range(0,height):
    for col in range(0,width):
        print str(row)+","+str(col)+":"+str(array.item(row,col))

Ale jeśli chcesz uzyskać współrzędną piksela, NumPyArray nie może ci pomóc. Możesz przekonwertować raster na punkt za pomocą narzędzia RasterToPoint, a następnie możesz uzyskać współrzędne według pola Kształt.

Najprostszym sposobem na wyprowadzenie współrzędnych i wartości komórek do pliku tekstowego w ArcGIS 10 jest funkcja przykładowa , bez potrzeby kodu, a zwłaszcza bez potrzeby zapętlania każdej komórki. W kalkulatorze rastrowym ArcGIS <= 9,3x był tak prosty, jak to, outfile.csv = sample(someraster)że generowałby plik tekstowy zawierający wszystkie wartości (nie zerowe) komórek i współrzędne (w formacie z, x, y). W ArcGIS 10 wygląda na to, że argument „in_location_data” jest teraz obowiązkowy, więc musisz użyć składni Sample(someraster, someraster, outcsvfile).

Edit: Można również określić wiele rastrów: Sample([someraster, anotherraster, etc], someraster, outcsvfile). Czy to zadziała na 8 miliardach komórek, nie mam pojęcia ...

Edycja: Uwaga, nie testowałem tego w ArcGIS 10, ale od lat korzystałem z przykładowej funkcji w <= 9,3 (i stacji roboczej).

Edycja: Testowałem już w ArcGIS 10 i nie będzie on wysyłany do pliku tekstowego. Narzędzie automatycznie zmienia rozszerzenie pliku na „.dbf”. Jednak ... następujący kod python działa, ponieważ instrukcje ArcMIS 10 są nadal obsługiwane przez instrukcje algebry map SOMA i MOMA :

import arcgisscripting
gp=arcgisscripting.create()
gp.multioutputmapalgebra(r'%s=sample(%s)' % (outputcsv,inputraster))

Jednym ze sposobów jest użycie narzędzia Raster_To_Point , a następnie narzędzia Add_XY_Coordinates . Otrzymasz plik kształtu, w którym każdy wiersz w tabeli atrybutów reprezentuje piksel z twojego rastra z kolumnami dla X_Coord , Y_Coord i Cell_Value . Następnie możesz zapętlić tę tabelę za pomocą kursora (lub wyeksportować ją do czegoś takiego jak Excel, jeśli wolisz).

Jeśli masz tylko jeden raster do przetworzenia, prawdopodobnie nie warto pisać skryptów - po prostu skorzystaj z narzędzi ArcToolbox. Jeśli musisz to zrobić dla wielu rastrów, możesz spróbować czegoś takiego:

[ Uwaga: nie mam ArcGIS 10 i nie znam ArcPy, więc jest to bardzo ogólny zarys. Jest nieprzetestowany i prawie na pewno będzie wymagał modyfikacji, aby go uruchomić.]

import arcpy, os
from arcpy import env

# User input
ras_fold = r'path/to/my/data'           # The folder containing the rasters
out_fold = r'path/to/output/shapefiles' # The folder in which to create the shapefiles

# Set the workspace
env.workspace = ras_fold

# Get a list of raster datasets in the raster folder
raster_list = arcpy.ListRasters("*", "All")

# Loop over the rasters
for raster in raster_list:
    # Get the name of the raster dataset without the file extension
    dataset_name = os.path.splitext(raster)[0]

    # Build a path for the output shapefile
    shp_path = os.path.join(out_fold, '%s.shp' % dataset_name)

    # Convert the raster to a point shapefile
    arcpy.RasterToPoint_conversion(raster, shp_path, "VALUE")

    # Add columns to the shapefile containing the X and Y co-ordinates
    arcpy.AddXY_management(shp_path)

Następnie można zapętlić tabele atrybutów shapefile za pomocą Search Cursor lub (być może prostszego) za pomocą dbfpy . Umożliwi to odczyt danych z twojego rastra (teraz zapisanego w tabeli shapefile .dbf) do zmiennych python.

from dbfpy import dbf

# Path to shapefile .dbf
dbf_path = r'path\to\my\dbf_file.dbf'

# Open the dbf file
db = dbf.Dbf(dbf_path)

# Loop over the records
for rec in db:
    cell_no = rec['POINTID'] # Numbered from top left, running left to right along each row
    cell_x = rec['POINT_X']
    cell_y = rec['POINT_Y']
    cell_val = rec['GRID_CODE']

    # Print values
    print cell_no, cell_x, cell_y, cell_val

Może mógłbyś utworzyć plik światowy dla rastra, ukryć raster w tablicy numpy. następnie, jeśli zapętlisz się nad tablicą, otrzymasz wartości komórek, a jeśli stopniowo zaktualizujesz x, y z pliku świata, będziesz również mieć współrzędne dla każdej wartości komórki. mam nadzieję, że się przyda.

Kod Marcina działał dobrze, z wyjątkiem tego, że problem z funkcjami rasCentrX i rasCentrY powodował, że współrzędne wyjściowe pojawiały się w innej rozdzielczości (jak zaobserwowała Grazia). Moją poprawką była zmiana

coordX = rasXmin + (0.5*rasMeanCellWidth + rasWidth)

do

coordX = rasXmin + ((0.5 + rasWidth) * rasMeanCellWidth)

i

  coordY = rasYmin + (0.5*rasMeanCellHeight + rasHeight)

do

  coordY = rasYmin + ((0.5 + rasHeight) * rasMeanCellHeight)

Użyłem tego kodu do konwersji siatki ESRI do pliku CSV. Osiągnięto to poprzez usunięcie odwołania do inRaster2, a następnie użycie csv.writer do wyprowadzenia współrzędnych i wartości:

out = csv.writer(open(outputCSV,"wb"), delimiter=',', quoting=csv.QUOTE_NONNUMERIC)
out.writerow(['X','Y','Value'])
(height, width, dim )=inRasterFullArray.shape
for row in range(0,height):
    for col in range(0,width):
        out.writerow(inRasterFullArray[row,col])

Nie znalazłem też potrzeby transpozycji

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray.T))

więc przekonwertowałem to na

inRasterFullArray = numpy.dstack((inRasterCoordinates, inRasterArray))

Brzydki, ale bardzo skuteczny:

1. Utwórz nowy element punktowy z 4 punktami poza narożnikami danego rastra. Upewnij się, że w tym samym układzie współrzędnych, co raster.
2. Dodaj podwójne pola „xcor” i „ycor”
3. Oblicz geometrię, aby uzyskać współrzędne tych pól
4. Analityk przestrzenny-> Interpolacja-> Trend -> Regresja liniowa
5. Ustawienia środowiska: przyciągnij raster i rozmiar komórki do tego samego, co dany raster
6. Wykonaj osobno dla „xcor” i „ycor”
7. Wychodzą miotacze ze współrzędnymi jako wartościami komórek, używanymi jako dane wejściowe dla skryptów.

Proste rozwiązanie z wykorzystaniem pakietów Pythona typu open source:

import fiona
import rasterio
from pprint import pprint


def raster_point_coords(raster, points):

    # initialize dict to hold data
    pt_data = {}

    with fiona.open(points, 'r') as src:
        for feature in src:
            # create dict entry for each feature
            pt_data[feature['id']] = feature

    with rasterio.open(raster, 'r') as src:
        # read raster into numpy array
        arr = src.read()
        # rasterio always reads into 3d array, this is 2d, so reshape
        arr = arr.reshape(arr.shape[1], arr.shape[2])
        # get affine, i.e. data needed to work between 'image' and 'raster' coords
        a = src.affine

    for key, val in pt_data.items():
        # get coordinates
        x, y = val['geometry']['coordinates'][0], val['geometry']['coordinates'][1]
        # use affine to convert to row, column
        col, row = ~a * (x, y)
        # remember numpy array is indexed array[row, column] ie. y, x
        val['raster_value'] = arr[int(row), int(col)]

    pprint(pt_data) 

if __name__ == '__main__':
    # my Landsat raster
    ras = '/data01/images/sandbox/LT05_040028_B1.tif'
    # my shapefile with two points which overlap raster area
    pts = '/data01/images/sandbox/points.shp'
    # call function
    raster_point_coords(ras, pts)

Fiona jest przydatna, ponieważ możesz otworzyć plik kształtu, iterować funkcje i (jak mam) dołączyć je do dictobiektu. Rzeczywiście, featuresama Fiona również jest podobna dict, więc łatwo jest uzyskać dostęp do właściwości. Gdyby moje punkty miały atrybuty, pojawiałyby się w tym dykcie wraz ze współrzędnymi, identyfikatorem itp.

Rasterio jest przydatny, ponieważ jest łatwy do odczytu w rastrze jako tablica liczb liczbowych, lekki i szybki typ danych. Mamy również dostęp do jednej dictz właściwości rastra, w tym do affinewszystkich danych potrzebnych do przekształcenia współrzędnych x, y rastra na współrzędne wiersza tablic i kolumn. Zobacz doskonałe wyjaśnienie @ perrygeo tutaj .

W efekcie powstaje pt_datatyp, dictktóry zawiera dane dla każdego punktu i jest wyodrębniany raster_value. Gdybyśmy chcieli, moglibyśmy również łatwo przepisać plik kształtu z wyodrębnionymi danymi.