Преобразование долготы, широты в значения пикселей с использованием GDAL - PullRequest
Новая
       75

Преобразование долготы, широты в значения пикселей с использованием GDAL

0 голосов
/ 12 октября 2018

У меня есть спутниковое изображение GeoTIFF и соответствующий файл OSM только с автомагистралями.Я хочу преобразовать значение широты долготы в файле OSM в пиксели и хочу выделить шоссе на спутниковом изображении.

Я пробовал несколько методов, которые описаны в StackExchange.Но я получаю отрицательное и одинаковое значение пикселя для каждого значения долготы и широты.Может кто-нибудь объяснить, что мне не хватает?

Вот информация об изображении, которое я собрал, используя приложение OTB.Image Description

Вот код, который я использую. 

from osgeo import gdal, osr
import numpy as np
import xml.etree.ElementTree as xml

src_filename = 'image.tif'
dst_filename = 'foo.tiff'


def readLongLat(path):
    lonlatList = []
    latlongtuple = ()
    root = xml.parse(path).getroot()
    for i in root:
        if i.tag == "node":
            latlong = []
            lat = float(i.attrib["lat"])
            long = float(i.attrib["lon"])
            latlong.append(lat)
            latlong.append(long)
            lonlatList.append(latlong)

    return lonlatList

# Opens source dataset
src_ds = gdal.Open(src_filename)
format = "GTiff"
driver = gdal.GetDriverByName(format)

# Open destination dataset
dst_ds = driver.CreateCopy(dst_filename, src_ds, 0)

# Get raster projection
epsg = 4269         # http://spatialreference.org/ref/sr-org/lambert_conformal_conic_2sp/
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)

# Make WGS84 lon lat coordinate system
world_sr = osr.SpatialReference()
world_sr.SetWellKnownGeogCS('WGS84')

transform = src_ds.GetGeoTransform()
gt = [transform[0],transform[1],0,transform[3],0,-transform[5]]

#Reading the osm file
lonlat = readLongLat("highways.osm")
# Transform lon lats into XY
coord_transform = osr.CoordinateTransformation(world_sr, srs)
newpoints = coord_transform.TransformPoints(lonlat) # list of XYZ tuples

# Make Inverse Geotransform  (try:except due to gdal version differences)
try:
    success, inverse_gt = gdal.InvGeoTransform(gt)
except:
    inverse_gt = gdal.InvGeoTransform(gt)

# [Note 1] Set pixel values
marker_array_r = np.array([[255]], dtype=np.uint8)
marker_array_g = np.array([[0]], dtype=np.uint8)
marker_array_b = np.array([[0]], dtype=np.uint8)
for x,y,z in newpoints:
    pix_x = int(inverse_gt[0] + inverse_gt[1] * x + inverse_gt[2] * y)
    pix_y = int(inverse_gt[3] + inverse_gt[4] * x + inverse_gt[5] * y)
    dst_ds.GetRasterBand(1).WriteArray(marker_array_r, pix_x, pix_y)
    dst_ds.GetRasterBand(2).WriteArray(marker_array_g, pix_x, pix_y)
    dst_ds.GetRasterBand(3).WriteArray(marker_array_b, pix_x, pix_y)

# Close files
dst_ds = None
src_ds = None

Ответы [ 2 ]

0 голосов
/ 15 октября 2018

Мне удалось это сделать с помощью библиотеки geoio. 

import geoio  

img = geoio.GeoImage(src_filename)
pix_x, pix_y = img.proj_to_raster(lon,lat)
0 голосов
/ 12 октября 2018

Недавно я попробовал использовать модуль xarray.Я думаю о xarray как о гибриде между pandas и numpy, который позволяет вам хранить информацию в виде массива, но обращаться к ней, используя просто .sel запросы.Документы здесь .

ОБНОВЛЕНИЕ: Похоже, что rasterio и xarray должны быть установлены для работы описанного ниже метода.См. ссылка .

Это гораздо более простой способ преобразования файла GeoTiff в удобный для пользователя массив.См. Мой пример ниже: 

import xarray as xr

ds = xr.open_rasterio("/path/to/image.tif")

# Insert your lat/lon/band below to extract corresponding pixel value
ds.sel(band=2, lat=19.9, lon=39.5, method='nearest').values
>>> [10.3]

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

Примечание: очевидно, нужно быть внимательным, чтобы убедиться, что ваши пары широта / долгота находятся в той же системе координат, что и файл GeoTiff, но я думаю, что вы все равно это обрабатываете.

...