Back to DrivenData | Blog

OpenAI carribean challenge lat / long to pixels

There is no topic on the forum for this competition. Can someone from please create it?
As usual in such competitions, at first you need to convert the geotiff pixels to latitude / longitude.

I remember successfully doing so with code like this:

def polygon_to_tuples(polygon):
    x, y = polygon.exterior.coords.xy
    x = list(x)
    y = list(y)
    return list(zip(x, y))

def coordinates_to_pixels(x, y, transform):
    """
    The Affine object is a named tuple with elements a, b, c, d, e, f
    corresponding to the elements in the matrix equation below,
    in which a pixel’s image coordinates are x, y and its world coordinates are x', y'.:

    | x' |   | a b c | | x |
    | y' | = | d e f | | y |
    | 1  |   | 0 0 1 | | 1 |
    """    
    rev = ~transform
    return rev * (x, y)

polygon = df.loc[100].geometry
tuples = polygon_to_tuples(polygon)
pixels = [coordinates_to_pixels(tup[0], tup[1], transform) for tup in tuples]

The coordinates of individual polygons seem to be more or less in line with house coordinates on Google Maps, but when I try to convert them to geo-tiff pixels, I get strange values.

Has anyone faced this issue?
I also tried using geoio with similar results.

1 Like

@aveysov when I use that same method to convert, I was interpreting them as meter locations, not pixel indices into the tiff image itself. But, this is my first exposure to working with TIFF files in this way, so I’m not sure that I’m interpreting it correctly. For reference, though, I have been using rasterio and when I convert the GeoJson coordinates using the same method that you do, I get values that seem to line up with the dataset.bounds property that is referenced here: https://rasterio.readthedocs.io/en/stable/quickstart.html#dataset-georeferencing.

That being said, I still get errors when I try to use the converted coordinates to mask the TIFF file and pull out a given ROI (although I think I might just have some x, y flipping that I need to do). Here’s a sample of the code I’ve been working with:

import json
import rasterio
from rasterio.mask import mask


fpath_tiff = '/data/caribbean_challenge/stac/colombia/borde_rural/borde_rural_ortho-cog.tif'
fpath_geojson = '/data/caribbean_challenge/stac/colombia/borde_rural/train-borde_rural.geojson'

with open(fpath_geojson) as geojson:
    geojson = json.load(geojson)

polygons = []
for feature in geojson['features']:
    polygons.append(feature['geometry'])
polygon0 = polygons[0]

def transform_coordinates(coordinates_lists, affine_transform):
    """Transform the coordinates with the provided `affine_transform`
    
    :param coordinates: list of lists of coordinates
    :type coordinates: list[list[int]]
    :param affine_transform: affine transformation to apply
    :type: affine.Affine
    """
    
    transformed_coordinates_lists = []
    for coordinates_list in coordinates_lists:
        transformed_coordinates_list = []
        for coordinate in coordinates_list:
            coordinate = tuple(coordinate)
            transformed_coordinate = affine_transform * coordinate
            transformed_coordinates_list.append(transformed_coordinate)
        transformed_coordinates_lists.append(transformed_coordinates_list)
    
    return [transformed_coordinates_lists]
    
unread_tiff = rasterio.open(fpath_tiff)
polygon0['coordinates'] = transform_coordinates(polygon0['coordinates'], unread_tiff.transform)

with rasterio.open(fpath_tiff) as tiff:
    out_image, out_transform = mask(tiff, [polygon0], crop=True)

The unread_tiff.bounds gives BoundingBox(left=592371.56518, bottom=501705.04172000004, right=593606.62878, top=503768.46364000003), which seems to be roughly in line with the converted lat / long values that I’m getting.

@aveysov not sure at the exact issue you’re having, but I was able to solve my issue by making sure that the GeoJSON lat / long coordinates were projected to the same CRS as the TIFF. The specific TIFF I had been playing with was in EPSG:32618, so I just used the pyproj to project the coordinates to that CRS rather than using the transform provided by the loaded TIFF. So, this seemed to work for me:

# Put together a comprehensive code example
fpath_tiff = '/data/caribbean_challenge/stac/colombia/borde_rural/borde_rural_ortho-cog.tif'
fpath_geojson = '/data/caribbean_challenge/stac/colombia/borde_rural/train-borde_rural.geojson'

import json

import rasterio
from pyproj import Proj
from rasterio.mask import mask
from rasterio.plot import show

%matplotlib inline

with open(fpath_geojson) as geojson:
    geojson = json.load(geojson)
    
polygons = []
for feature in geojson['features']:
    polygons.append(feature['geometry'])
polygon0 = polygons[0]

def transform_coordinates(coordinates_lists, transform):
    """Transform the coordinates with the provided `affine_transform`
    
    :param coordinates_lists: list of lists of coordinates
    :type coordinates_lists: list[list[int]]
    :param transform: transformation to apply
    :type: pyproj.Proj
    """
    
    transformed_coordinates_lists = []
    for coordinates_list in coordinates_lists:
        transformed_coordinates_list = []
        for coordinate in coordinates_list:
            coordinate = tuple(coordinate)
            transformed_coordinate = list(transform(coordinate[0], coordinate[1]))
            transformed_coordinates_list.append(transformed_coordinate)
        transformed_coordinates_lists.append(transformed_coordinates_list)
    
    return transformed_coordinates_lists
    
unread_tiff = rasterio.open(fpath_tiff)
proj = Proj(init='epsg:32618')
polygon0['coordinates'] = transform_coordinates(polygon0['coordinates'], proj)

with rasterio.open(fpath_tiff) as tiff:
    out_image, out_transform = mask(tiff, [polygon0], crop=True)

show(out_image)

Resulting image:

2 Likes

Hi @sallamander , looks like this solves my problem.
I will be testing it tomorrow, but that looks like the missing piece of the puzzle

1 Like

Changed your example a bit to rely as less as possible on the geo-stack libraries, which are poorly documented in my opinion.

So, I hope this helps someone (also it is much easier to iterate over the whole dataset in such fashion, using plain vanilla data-structures like lists and pandas dataframes):

import rasterio
import pandas as pd
from pyproj import Proj
import geopandas as gpd
from pathlib import Path

curr_path = Path('')
data_path = curr_path / '../data'
stac_path = data_path / 'stac'

sub_df = pd.read_csv('../data/submission_format.csv')
trn_df = pd.read_csv('../data/train_labels.csv')
met_df = pd.read_csv('../data/metadata.csv')

tiff_paths = list(met_df.image.values)
train_geo_paths = list(met_df.train.values)
test_geo_paths = list(met_df.test.values)


def polygon_to_tuples(polygon):
    x, y = polygon.exterior.coords.xy
    x = list(x)
    y = list(y)
    return list(zip(x, y))


def transform_coordinates(coordinate_lists,
                          transform):
    transformed_coordinates_lists = []
    for coordinates_list in coordinate_lists:
        transformed_coordinates_list = []
        for coordinates in coordinates_list:
            transformed_coordinate = tuple(transform(coordinates[0],
                                                     coordinates[1]))
            transformed_coordinates_list.append(transformed_coordinate)
        transformed_coordinates_lists.append(transformed_coordinates_list)
    return transformed_coordinates_lists


def get_train_df(train_geo_paths,
                 tiff_paths,
                 data_path):
    dfs = []

    for train_path, tiff_path in zip(train_geo_paths,
                                     tiff_paths):
        train_path = str(data_path / train_path)
        tiff_path = str(data_path / tiff_path)
        alias = train_path.split('/')[-2]
        df = gpd.read_file(train_path)
        df['train_path'] = train_path
        df['tiff_path'] = tiff_path
        df['alias'] = alias
        with rasterio.open(tiff_path) as tiff:
            epsg = tiff.crs.to_epsg()
            df['epsg'] = epsg
        dfs.append(df)

    df = pd.concat(dfs)
    
    return df


train_df = get_train_df(train_geo_paths,
                        tiff_paths,
                        data_path)


sample = train_df.iloc[0]
polygon = sample.geometry
tuples = polygon_to_tuples(polygon)
proj = Proj(init='epsg:{}'.format(sample.epsg))
tiff_path = sample.tiff_path
coordinates = transform_coordinates([tuples], proj)[0]

with rasterio.open(tiff_path) as tiff:
    pixels = [tiff.index(*coord) for coord in coordinates]
1 Like

Oh, there are multi-polygons there as well
Maybe use x, y = polygon.convex_hull.exterior.coords.xy for them

Hi @sallamander, thanks for your code. It has helped me a lot.

I’m still having an error for a few points on both train and test. The message I get is: “ValueError: Input shapes do not overlap raster.” and I was wondering if you have solved it.

It should be related to the CRS as well, but I have made sure to use the same that the TIFF is using (epsg:32616).
These are some of the IDs where the issue happens: 7a305846, 7a277410, 7a39f572.

Thanks a lot!

maybe @aveysov has seen this issue as well?

@javiergalvis glad to hear it helped! After sending that along I actually ended up pivoting to use GeoPandas, which appears to handle all the conversions for you. When using it I didn’t run into any errors across the training or test sets when pulling the roof pixels out of the TIFFs themselves. I’m pulling this code out of a class so it might not work exactly as is (even after you fill in the filepaths), but should capture the jist of it:

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.mask import mask

fpath_tiff = # put filepath of tiff here 
fpath_geojson = # put filepath of geojson here

df_roof_geometries = gpd.read_file(fpath_geojson)
roof_geometries = (
    df_roof_geometries[['id', 'projected_geometry']].values
)

with rasterio.open(fpath_tiff) as tiff:
            for roof_id, projected_geometry in roof_geometries:
                roof_image, _ = mask(
                    tiff, [projected_geometry], crop=True, pad=True,
                    filled=False, pad_width=padding
                )
                roof_image = np.transpose(roof_image, (1, 2, 0))
                roof_mask, _ = mask(
                    tiff, [projected_geometry], crop=True, pad=True,
                    filled=True, pad_width=padding
                )
                roof_mask = np.transpose(roof_mask, (1, 2, 0))
                roof_mask = roof_mask[..., 0] > 0
                
                # do something with roof_id, roof_image, and / or roof_mask

Let me know if you have any more questions about this code or I can help further with any debugging. Hopefully this helps!

Thanks a lot for replying!
how are you projecting your geometries? I realized my issue happens only with multipolygons, polygons are being masked without any problem.

Oh, whoops I forgot that bit, sorry about that! Here’s the relevant code snippet:

        with rasterio.open(fpath_tiff) as tiff:
            tiff_crs = tiff.crs.data
            df_roof_geometries['projected_geometry'] = (
                df_roof_geometries['geometry'].to_crs(tiff_crs)
            )

I think I did have the same original problem of having MultiPolygons not being projected and / or parsed properly, and then after switching to GeoPandas (mainly for usability reasons), it was resolved so I didn’t really dig any further. I’ve visualized a fair amount of the images that have been generated with that code and everything looks to be parsed correctly.

1 Like

amazing! thanks a lot!