Создание наноразмерной матрицы высот с помощью GDAL

Возможно, это немного странный вопрос, но позвольте мне дать вам краткое объяснение истории вопроса перед моими фактическими вопросами:

Атомно-силовая микроскопия (АСМ) - это метод, который, вкратце (и насколько мне известно), позволяет исследователям сканировать области на микро- и наноразмерном уровне. Он работает путем «сканирования» области с помощью своего рода зонда. Мне сложно объяснить еще кое-что, поскольку я не понимаю этого по-настоящему. Что я действительно знаю, и что вызвало мое любопытство, так это то, что в результате получается «сетка» значений «высоты» (матрица, скажем, значений 512x512, описывающая высоту зонда в этой точке).

Затем я подумал: ну, кроме масштаба, это на самом деле цифровая модель рельефа! И это означает, что если мне удастся создать файл DEM, как это понимают ГИС-инструменты, я смогу применить к нему ГИС-анализ!

Как бы то ни было, моя другая работа в лаборатории, где есть АСМ-машина, и использует ее в одном из своих проектов. Я получил от нее несколько файлов сканирования, и мне удалось, используя Python (struct и numpy), проанализировать эти двоичные файлы, и теперь у меня есть пустой массив размером 512x512, заполненный значениями int16.

То, что я планирую дальше, и в чем мне нужна некоторая помощь, это «отображение на правильную матрицу высот». У меня есть некоторые знания о DEMS, но когда дело доходит до их фактического поколения, я совершенно новый.

Я думаю, что мне нужно как-то привязать мои данные, и для этого мне нужна пользовательская (плоская) система координат. Я предполагаю, что моя система координат будет использовать микро- или нанометры в качестве единиц измерения. Тогда нужно просто определить размер области, отсканированной с помощью AFM (я полагаю, это где-то в двоичном файле, предположим, что это известно).

обновление : у меня также есть несколько сканов в разных разрешениях, но одной и той же области. Например, у меня есть эта информация о двух сканированиях:

увеличенное изображение:

Scan Size: 51443.5 nm
X Offset: 0 nm
Y Offset: 0 nm

уменьшенное (детальное) изображение:

Scan Size: 5907.44 nm
X Offset: 8776.47 nm
Y Offset: 1486.78 nm

Я имею в виду, что моя пользовательская система координат должна иметь начало в 0,0, а для большего изображения я могу назначить пикселю 0,0 значение координаты (0,0), а пикселю 512,512 значение координаты (51443.5, 51443.5 ) (Думаю, вы получите картину для других необходимых точек).

Затем большее изображение будет отображать пиксели (0,0) в (8776,47, 1486,78) и (512,512) в (8776,47 + 5907,44, 1486,78 + 5907,44)

Первый вопрос : как мне создать proj4 def для такой системы координат? То есть: как я могу назначить эти «координаты реального мира» моей пользовательской системе координат (или, если я последую совету Уубера и использую локальную систему координат и лгу о единицах (т.е. рассматриваю мои нанометры как километры)

Затем я должен перенести свой массив 2-мерного массива в формат файлов с географической привязкой DEM. Я думал об использовании GDAL (или, скорее, привязки Python).

Второй вопрос : как мне создать DEM с географической привязкой из «произвольных» данных, таких как мои? Желательно в Python и с использованием библиотек с открытым исходным кодом.

Остальное должно быть достаточно легким делом, просто вопрос использования правильных инструментов анализа. Проблема в том, что эта задача обусловлена моим собственным любопытством, поэтому я не совсем уверен, что мне на самом деле нужно делать с наноразмерной ЦМР. Это напрашивается на

Третий вопрос : что делать с наноразмерной ЦМР? Какого рода анализ может быть выполнен, каковы подходящие инструменты для анализа матрицы высот и, наконец: возможно ли сделать карту с отмывками и контурными линиями из этих данных? :)

Я приветствую все предложения и указания, но имейте в виду, что я ищу бесплатные альтернативы, так как это строго основанный на хобби проект без бюджета или финансирования (и у меня нет доступа к любым лицензированным ГИС-приложениям). Кроме того, я знаю, что Bruker, компания, которая продает эти машины AFM, поставляет некоторое программное обеспечение, но использовать его было бы неинтересно.

python coordinate-system dem atlefren
источник

Весело и интересно! Можете ли вы опубликовать пример данных? Требуется ли нанометровая шкала на проекции? Думая, может быть, легче масштабироваться, хотя это немного "обманывает". Между прочим, я думаю, что вы можете пройти долгий путь с помощью GDAL / ogr, хотя проблему с проекцией все равно придется решать. gdal.org/gdal_grid.html

alexanno

Благодарность! Думаю, это больше комментарий, чем ответ. Что касается нанометрового масштаба, я бы сказал, что все, что работает, замечательно, но истинная наноразмерная проекция была бы самой крутой. Когда дело доходит до выборки данных, мне нужно проверить, есть ли какие-то ограничения, но в основном это 2-мерная матрица значений int16.

atlefren

Зачем вам нужны параметры proj4? Вы не будете передавать эти данные в другую систему координат (особенно географическую). Имхо, вы должны быть в состоянии сделать весь свой анализ без какой-либо системы координат. Что вы понимаете под DEM? Существует несколько типов, таких как триангулированные поверхности (например, триангуляция Делона) или растровые карты (у вас это уже есть). Это, конечно, сильно зависит от программного обеспечения для анализа. Конечно, вы можете создавать другие карты, если они необходимы в качестве вывода, чтобы понять зонд. Вы можете взглянуть на code.google.com/p/mtex для анализа зерна.

ошибка

Причина, по которой я (думаю) мне нужен CRS, заключается в следующем: если я просто создаю, скажем, GeoTIFF без назначения CRS, единицей измерения будут пиксели. Что если я хочу измерить расстояния? И что, если у меня есть два AFM-сканирования и я знаю, как они связаны друг с другом (с точки зрения масштаба и смещения от некоторой точки). Назначение CRS облегчило бы просмотр нескольких сканов одновременно

atlefren

Я обычно справляюсь с такими данными, устанавливая локальную систему координат (с началом координат, совпадающим с источником изображения) и «лгая» относительно единиц измерения. Например, вы могли бы оговорить, что единицы измерения - это километры, когда они действительно являются нанометрами, что облегчает умственное преобразование туда и обратно. Конечно, вы не будете заниматься перепроектированием, так что это не проблема. Настройка этой системы координат такая же, как и географическая привязка любой матрицы высот; это может быть так же просто, как создать файл без поворота.

whuber

Ну, похоже, я решил проблемы 1 и 2, по крайней мере. Полный исходный код на github , но некоторые пояснения здесь:

Для пользовательского CRS я решил (согласно предложению Whubers) «обмануть» и использовать счетчики в качестве единицы измерения. Я нашел «локальные crs» на apatialreference.org ( SR-ORG: 6707 ):

LOCAL_CS["Non-Earth (Meter)",
    LOCAL_DATUM["Local Datum",0],
    UNIT["Meter",1.0],
    AXIS["X",NORTH],
    AXIS["Y",EAST]]

Используя Python и GDAL, это довольно легко прочитать:

def get_coordsys():
    #load our custom crs
    prj_text = open("coordsys.wkt", 'r').read()
    srs = osr.SpatialReference()
    if srs.ImportFromWkt(prj_text):
        raise ValueError("Error importing PRJ information" )

    return srs

Кроме того, написание DEM с помощью GDAL было на самом деле довольно простым (я закончил с однополосным геотиффом). Строка parser.read_layer (0) возвращает мою ранее описанную матрицу 512x512.

def create_dem(afmfile, outfile):

    #parse the afm data
    parser = AFMParser(afmfile)

    #flip to account for the fact that the matrix is top-left, but gdal is bottom-left
    data = flipud(parser.read_layer(0))

    driver = gdal.GetDriverByName("GTiff")
    dst_ds = driver.Create(
        outfile,
        data.shape[1],
        data.shape[0],
        1 ,
        gdal.GDT_Int16 ,
    )

    dst_ds.SetGeoTransform(get_transform(parser, data))
    dst_ds.SetProjection(get_coordsys().ExportToWkt())
    dst_ds.GetRasterBand(1).WriteArray(data)
    dst_ds = None

Самая трикотажная часть заключалась в том, чтобы выяснить, как правильно «привязать» мой файл, в итоге я использовал SetGeoTransform , получив параметры следующим образом:

def get_transform(parser, data):
    #calculate dims
    scan_size, x_offset, y_offset = parser.get_size()
    x_size = data.shape[0]
    y_size = data.shape[1]
    x_res = scan_size / x_size
    y_res = scan_size / y_size

    #set the transform (offsets doesn't seem to work the way I think)
    #top left x, w-e pixel resolution, rotation, 0 if image is "north up", top left y, rotation, 0 if image is "north up", n-s pixel resolution
    return [0, x_res, 0, 0, 0, y_res]

Эта последняя часть, пожалуй, та, в которой я больше всего не уверен, что я действительно искал, это что-то вроде строки * gdal_transform -ullr *, но я не мог найти способ сделать это программно.

Я могу открыть свой GeoTIFF в Qgis и просмотреть его (и визуально сравнивая его с результатом программы Bruker, он выглядит правильно), но я на самом деле не ответил на мой вопрос 3; что делать с этими данными. Итак, я открыт для предложений!

atlefren
источник

Одна интересная задача может состоять в том, чтобы сравнить расстояния на ЦМР с расстояниями между точками земного шара, чтобы дать зрителям представление о том, насколько маленький наноразмер. Например, htwins.net/scale2

blah238

Создание наноразмерной матрицы высот с помощью GDAL

Ответы: