У меня есть двоичная сетка Arc / Info - в частности, растр накопления потока ArcGIS - и я хотел бы идентифицировать все ячейки, имеющие определенное значение (или диапазон значений). В конечном счете, мне нужен шейп-файл точек, представляющих эти ячейки.
Я могу использовать QGIS, чтобы открыть hdr.adf и получить этот результат, рабочий процесс:
- QGIS> Растровое меню> Растровый калькулятор (отметьте все точки целевым значением)
- QGIS> Растровое меню> Полигонизация
- QGIS> Векторное меню> Подменю геометрии> Центроиды многоугольников
- Отредактируйте центроиды, чтобы удалить ненужные полицентроиды (те = 0)
Этот подход «выполняет свою работу», но он мне не нравится, потому что он создает 2 файла, которые мне нужно удалить, а затем я должен удалить ненужные записи из шейп-файла центроидов (то есть те, которые = 0).
Существующий вопрос подходит к этой теме, но это специально для ArcGIS / ArcPy, и я хотел бы остаться в пространстве FOSS.
Есть ли у кого-нибудь существующий рецепт / сценарий GDAL / Python, который запрашивает значения ячеек растра, и когда целевое значение - или значение в целевом диапазоне - найдено, запись добавляется в шейп-файл? Это не только позволит избежать взаимодействия с пользовательским интерфейсом, но и даст чистый результат за один проход.
Я сделал снимок, работая против одной из презентаций Криса Гаррарда , но растровых работ нет в моей рулевой рубке, и я не хочу загромождать вопрос своим слабым кодом.
Если кто-то захочет поиграть с точным набором данных, я поставлю его здесь как .zip .
[Редактировать заметки] Оставив это для потомков. Смотрите обмен комментариями с om_henners. В основном значения x / y (строка / столбец) были перевернуты. Оригинальный ответ содержал эту строку:
(y_index, x_index) = np.nonzero(a == 1000)
перевернутый, вот так:
(x_index, y_index) = np.nonzero(a == 1000)
Когда я впервые столкнулся с проблемой, показанной на скриншоте, я подумал, правильно ли я реализовал геометрию, и я экспериментировал, переключая значения координат x / y в этой строке:
point.SetPoint(0, x, y)
..в качестве..
point.SetPoint(0, y, x)
Однако это не сработало. И я не думал пытаться перевернуть значения в выражении Numpy от om_henners, ошибочно полагая, что переворачивание их в любой строке было эквивалентно. Я считаю , что реальная проблема относится к x_size
и y_size
значениям, соответственно 30
и -30
, которые применяются , когда строки и столбцов индексы используются для вычисления координат точек для ячеек.
[Исходное редактирование]
@om_henners, я пробую ваше решение в сочетании с парой получателей для создания шейп-файлов точек с помощью ogr ( invisibleroads.com , Крис Гаррард ), но у меня возникла проблема, когда точки отображаются так, как если бы они отражались через проход до 315/135 градусов.
Светло-голубые точки : созданные моим подходом QGIS , выше
Фиолетовые точки : созданные с помощью кода Py GDAL / OGR , ниже
[Решаемые]
Этот код Python реализует полное решение, предложенное @om_henners. Я проверил это, и это работает. Спасибо чувак!
from osgeo import gdal
import numpy as np
import osgeo.ogr
import osgeo.osr
path = "D:/GIS/greeneCty/Greene_DEM/GreeneDEM30m/flowacc_gree/hdr.adf"
print "\nOpening: " + path + "\n"
r = gdal.Open(path)
band = r.GetRasterBand(1)
(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = r.GetGeoTransform()
a = band.ReadAsArray().astype(np.float)
# This evaluation makes x/y arrays for all cell values in a range.
# I knew how many points I should get for ==1000 and wanted to test it.
(y_index, x_index) = np.nonzero((a > 999) & (a < 1001))
# This evaluation makes x/y arrays for all cells having the fixed value, 1000.
#(y_index, x_index) = np.nonzero(a == 1000)
# DEBUG: take a look at the arrays..
#print repr((y_index, x_index))
# Init the shapefile stuff..
srs = osgeo.osr.SpatialReference()
#srs.ImportFromProj4('+proj=utm +zone=15 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')
srs.ImportFromWkt(r.GetProjection())
driver = osgeo.ogr.GetDriverByName('ESRI Shapefile')
shapeData = driver.CreateDataSource('D:/GIS/01_tutorials/flow_acc/ogr_pts.shp')
layer = shapeData.CreateLayer('ogr_pts', srs, osgeo.ogr.wkbPoint)
layerDefinition = layer.GetLayerDefn()
# Iterate over the Numpy points..
i = 0
for x_coord in x_index:
x = x_index[i] * x_size + upper_left_x + (x_size / 2) #add half the cell size
y = y_index[i] * y_size + upper_left_y + (y_size / 2) #to centre the point
# DEBUG: take a look at the coords..
#print "Coords: " + str(x) + ", " + str(y)
point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
point.SetPoint(0, x, y)
feature = osgeo.ogr.Feature(layerDefinition)
feature.SetGeometry(point)
feature.SetFID(i)
layer.CreateFeature(feature)
i += 1
shapeData.Destroy()
print "done! " + str(i) + " points found!"
srs.ImportFromWkt(r.GetProjection())
(вместо того, чтобы создавать проекцию из известной строки проекта).Ответы:
В GDAL вы можете импортировать растр в виде массива.
Затем, используя numpy, легко получить индексы массива, соответствующие булеву выражению:
Из растрового геотрансформации мы можем получить такую информацию, как координаты x и y в верхнем левом углу и размеры ячеек.
Верхняя левая ячейка соответствует
a[0, 0]
. Размер Y всегда будет отрицательным, поэтому, используя индексы x и y, вы можете рассчитать координаты каждой ячейки на основе индексов.Отсюда достаточно просто создать шейп-файл с помощью OGR. Пример кода приведен в этом вопросе о том, как создать новый набор данных с точечной информацией.
источник
(y_index, x_index) = np.nonzero(a > threshold)
.shp
--- только то, что не работа, и при этом это не было где-нибудь близко. Я не был шокирован, так как значение x исчисляется сотнями тысяч, а значение y - миллионами, поэтому меня это сильно смутило. Я не думал, чтобы попытаться перевернуть их обратно на выражение Numpy. Большое спасибо за вашу помощь, это круто. Именно то, что я хотел. :)Почему бы не использовать набор инструментов Sexante в QGIS? Это похоже на построитель моделей для ArcGIS. Таким образом, вы можете связывать операции и рассматривать их как одну операцию. Вы можете автоматизировать удаление промежуточных файлов и удаление ненужных записей, если я не ошибаюсь.
источник
Может быть полезно импортировать данные в postgis (с поддержкой растров) и использовать там функции. Этот урок может содержать элементы, которые вам нужны.
источник