Маскированный растр с помощью Shapefile с использованием растерио: как сохранить информацию шейп-файла, соответствующую замаскированному растровому пикселю - PullRequest
Новая
       46

Маскированный растр с помощью Shapefile с использованием растерио: как сохранить информацию шейп-файла, соответствующую замаскированному растровому пикселю

0 голосов
/ 23 октября 2019

У меня есть растр и шейп-файл, который содержит много полигонов (контуры здания). Каждый полигон в шейп-файле имеет ассоциированное значение (уровень урона от 0 до 3). Я создал замаскированный растр с помощью rasterio.mask.mask и преобразовал выходные данные в геоданный кадр (называемый 'd'), который содержит координаты широты / долготы центральной точки замаскированных растровых пикселей, связанного значения растрового пикселя (когерентностьзначение 'cc'). Тем не менее, я хочу, чтобы последний геодатафрейм содержал уровни повреждений и геометрию полигона из шейп-файла, используемого для создания маски. Как бы я это сделал? Мой код для создания маски приведен ниже (на основе https://gis.stackexchange.com/questions/260304/extract-raster-values-within-shapefile-with-pygeoprocessing-or-gdal):

import rasterio as rio
cc_path=(r'coherenceND_cut.tif')
cc_data=rio.open(cc_path)
out_image, out_transform=rio.mask.mask(cc_data, shapefile['geometry'], all_touched=True, nodata=999, crop=True)
out_meta = cc_data.meta

nodata=999
# extract the values of the masked array
extractdata = out_image[0,:,:]
# extract the row, columns of the valid values
import numpy as np
row, col = np.where(extractdata != nodata) 
cc = np.extract(extractdata != nodata, extractdata)

from rasterio import Affine # or from affine import Affine
T1 = out_transform * Affine.translation(0.5, 0.5) # reference the pixel centre
rc2xy = lambda r, c: (c, r) * T1  

d = gpd.GeoDataFrame({'col':col,'row':row,'cc':cc})
# coordinate transformation
d['x'] = d.apply(lambda row: rc2xy(row.row,row.col)[0], axis=1)
d['y'] = d.apply(lambda row: rc2xy(row.row,row.col)[1], axis=1)
# geometry
from shapely.geometry import Point
d['geometry'] =d.apply(lambda row: Point(row['x'], row['y']), axis=1)

Первые 2 строки кадра выходных данных (d) выглядят так: 

col row cc              x               y               
222 4   -0.0102040805   13.29583338772  42.77888886942

Но этохочу, чтобы это выглядело так: 

col row cc              x               y                    damage level    polygon
222 4   -0.0102040805   13.29583338772  42.77888886942     1             ((13.23572745800004 42.73806166700007 0, 13.23560465900005 42.73804756200008 0, 13.23560623600002 42.73808437900004 0, 13.23572871300007 42.73808739300006 0, 13.23572785700003 42.73806985400006 0, 13.23572745800004 42.73806166700007 0))
...