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.
2 Likes
@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]
3 Likes
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?
1 Like
@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!
4 Likes
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.
1 Like
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.
4 Likes
To add on to the other great comments in the thread, this is what I’ve been using to reproject lat/long ESPG:4326 to ESPG:<whatever_the_geotiff_is> then to row/col to index into the geotiff image:
geojson_data = geojson.load("path/to/geojson")
feature = geojson_data['features'][0]
with rasterio.open(geotiff_path) as ds:
src_proj = str(ds.crs)
t = ds.transform
src_crs = rasterio.crs.CRS.from_string('EPSG:4326')
tgt_crs = rasterio.crs.CRS.from_string(src_proj)
warped_geometry = rasterio.warp.transform_geom(src_crs, tgt_crs, feature['geometry'], precision=8)
Then once you have it reprojected into the geotiff’s projection you can get the row and column using:
row_coords, col_coords = (rasterio.transform.rowcol(t, x_coords, y_coords))
There has to be a one liner to do this in geopandas or something but I couldn’t figure it out so I just loop though all the features in the .geojson and update them. @sallamander solution is nice and clean. The part I got stuck on was I assumed the geotiff was in ESPG:4326 but it actually changes per image
stac/guatemala/mixco_1_and_ebenezer/mixco_1_and_ebenezer_ortho-cog.tif, EPSG:32616
stac/guatemala/mixco_3/mixco_3_ortho-cog.tif, EPSG:32616
stac/st_lucia/castries/castries_ortho-cog.tif, EPSG:32620
stac/st_lucia/dennery/dennery_ortho-cog.tif, EPSG:32620
stac/st_lucia/gros_islet/gros_islet_ortho-cog.tif, EPSG:32620
stac/colombia/borde_rural/borde_rural_ortho-cog.tif, EPSG:32618
stac/colombia/borde_soacha/borde_soacha_ortho-cog.tif, EPSG:32618
@aveysov How comes some of my pixel values are higher than 255. Which format is it under? Thanks 
@buildwithcynthia probably you use RGB with three channel. try to convert them in gray scale then check
Hey, just curious!!
@sallamander Why are you not visible on Leaderboard, did not proceed further (dumb question)?
