CODEWITHSUNDEEP
pythoncoordinatesgispixelgdal
user1703942
user1703942

Reputation: 327

Convert Longitude, latitude to Pixel values using GDAL

I have a satellite GeoTIFF Image and a corresponding OSM file with only the highways. I want to convert the longitude latitude value in the OSM file to the pixels and want to highlight highway on the satellite image.

I have tried several methods that are explained on StackExchange. But I get the negative and same pixel value for every longitude and latitude values. Could somebody explain, what am I missing?

Here is the information of the image that I have gathered using OTB application. Image Description

Here is the code that i am using.

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

Upvotes: 2

Views: 4441

Answers (2)

user1703942
user1703942

Reputation: 327

I was able to do that using the library geoio.

import geoio  

img = geoio.GeoImage(src_filename)
pix_x, pix_y = img.proj_to_raster(lon,lat)

Upvotes: 0

tda
tda

Reputation: 2133

Something I have tried recently is using the xarray module. I think of xarray as a hybrid between pandas and numpy that allows you to store information as an array but access it using simply .sel requests. Docs here.

UPDATE: Seems as if rasterio and xarray are required to be installed for the below method to work. See link.

It is a much simpler way of translating a GeoTiff file to a user-friendly array. See my example below:

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]

This does not answer your question directly, but may help you identify a different (and probably simpler) approach that I've recently switched to.

Note: obviously care needs to be taken to ensure that your lat/lon pairs are in the same coordinate system as the GeoTiff file, but I think you're handling that anyway.

Upvotes: 1

Related Questions