2

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
user1703942
  • 317
  • 3
  • 15

2 Answers2

1

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.

tda
  • 2,045
  • 1
  • 17
  • 42
  • Thank you for your reply. I am trying your method. But I got the error message that AttributeError: module 'xarray' has no attribute 'open_rasterio'. Although I have correctly installed Xarray. I am using Python 3.5.4. – user1703942 Oct 12 '18 at 12:31
  • What version of `xarray` do you have installed? Mine is 0.10.4. – tda Oct 12 '18 at 12:35
  • I had installed 0.10.9. But I reinstalled the version 0.10.4 but still have the same error. which Python version are you using? – user1703942 Oct 12 '18 at 12:40
  • v3.6.4. I didn't realise that would make a difference - damn. Only other thing is whether you have `rasterio` installed - could be the issue? Updated answer: `rasterio` is a requirement. – tda Oct 12 '18 at 13:06
  • That's interesting?! I'm afraid I can't replicate even when building a Python 3.5 environment from scratch... I get an error stating 'No module named 'rasterio' then I install `rasterio` and works fine. – tda Oct 12 '18 at 13:55
0

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)
user1703942
  • 317
  • 3
  • 15