Skip to content

Commit ede627d

Browse files
avanettenavanetten
committed
Update add_geo_coords.py
1 parent 93776e9 commit ede627d

File tree

1 file changed

+219
-12
lines changed

1 file changed

+219
-12
lines changed

‎core/add_geo_coords.py‎

Lines changed: 219 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,11 @@
55
66
@author: avanetten
77
8+
See:
89
https://github.com/SpaceNetChallenge/utilities/blob/spacenetV3/spacenetutilities/geoTools.py
9-
1010
https://gis.stackexchange.com/questions/127427/transforming-shapely-polygon-and-multipolygon-objects
11-
1211
"""
1312

14-
15-
import os
1613
import gdal
1714
import geopandas as gpd
1815
import shapely
@@ -23,15 +20,225 @@
2320
import pandas as pd
2421
import numpy as np
2522
import time
26-
#import json
23+
# import json
24+
# import os
25+
26+
27+
###############################################################################
28+
def geomPixel2geomGeo(shapely_geom, affineObject=[],
29+
input_raster='', gdal_geomTransform=[]):
30+
31+
"""
32+
Transform a shapely geometry from pixel to geospatial coordinates
33+
34+
Arguments
35+
---------
36+
shapely_geom : shapely.geometry
37+
Input geometry in pixel coordinates
38+
affineObject : affine Object ?
39+
Object that is used to transform geometries (Preferred method).
40+
This argument supercedes input_raster. Defaults to ``[]``.
41+
input_raster : str
42+
Path to raster to gather georectification information. If string is
43+
empty, use gdal_geomTransform. Defaults to ``''``.
44+
gdal_geomTransform : gdal transform?
45+
The geographic transform of the image. Defaults to ``[]``.
46+
47+
Returns
48+
-------
49+
geomTransform : shapely transform?
50+
Geometric transform.
51+
"""
52+
53+
if not affineObject:
54+
if input_raster != '':
55+
affineObject = rio.open(input_raster).transform
56+
else:
57+
affineObject = af.Affine.from_gdal(gdal_geomTransform)
58+
59+
if not gdal_geomTransform:
60+
gdal_geomTransform = gdal.Open(input_raster).GetGeoTransform()
61+
62+
geomTransform = shapely.affinity.affine_transform(shapely_geom,
63+
[affineObject.a,
64+
affineObject.b,
65+
affineObject.d,
66+
affineObject.e,
67+
affineObject.xoff,
68+
affineObject.yoff]
69+
)
70+
71+
return geomTransform
72+
73+
74+
###############################################################################
75+
def get_row_geo_coords(row, affineObject=[], gdal_geomTransform=[],
76+
inProj_str='epsg:4326', outProj_str='epsg:3857',
77+
verbose=False):
78+
79+
"""
80+
Get geo coords for SIMRSWN dataframe row containing pixel coords.
81+
82+
Arguments
83+
---------
84+
row : pandas dataframe
85+
Row in the pandas dataframe
86+
affineObject : affine Object ?
87+
Object that is used to transform geometries. Defaults to ``[]``.
88+
gdal_geomTransform : gdal transform?
89+
The geographic transform of the image. Defaults to ``[]``.
90+
inProj_str : str
91+
Projection of input image. Defaults to ``'epsg:4326'`` (WGS84)
92+
outProj_str : str
93+
Desired projection of bounding box coordinate.
94+
Defaults to ``'epsg:3857'`` (WGS84 Web Mercator for mapping).
95+
verbose : boolean
96+
Switch to print values to screen. Defaults to ``False``.
97+
98+
Returns
99+
-------
100+
out_arr, poly_geo : tuple
101+
out_arr is a list of geo coordinates computd from pixels coords:
102+
[lon0, lat0, lon1, lat1, x0_wmp, y0_wmp, x1_wmp, y1_wmp]
103+
poly_geo is the shapely polygon in geo coords of the box
104+
"""
105+
106+
# convert latlon to wmp
107+
inProj = pyproj.Proj(init=inProj_str) # (init='epsg:4326')
108+
outProj = pyproj.Proj(init=outProj_str) # init='epsg:3857')
109+
110+
raster_loc = row['Image_Path']
111+
x0, y0 = row['Xmin_Glob'], row['Ymin_Glob']
112+
x1, y1 = row['Xmax_Glob'], row['Ymax_Glob']
113+
if verbose:
114+
print("idx, x0, y0, x1, y1:", row.values[0], x0, y0, x1, y1)
115+
116+
poly = shapely.geometry.Polygon([(x0, y0), (x0, y1), (x1, y1), (x1, y0)])
117+
# transform
118+
t01 = time.time()
119+
poly_geo = geomPixel2geomGeo(poly, input_raster=raster_loc,
120+
affineObject=affineObject,
121+
gdal_geomTransform=gdal_geomTransform)
122+
# print ("geo_coords.coords:", list(poly_geo.exterior.coords))
123+
t02 = time.time()
124+
if verbose:
125+
print(" Time to compute transform:", t02-t01, "seconds")
126+
127+
# x-y bounding box is a (minx, miny, maxx, maxy) tuple.
128+
lon0, lat0, lon1, lat1 = poly_geo.bounds
129+
# wkt_latlon = poly_geo.wkt
130+
if verbose:
131+
print(" lon0, lat0, lon1, lat1:", lon0, lat0, lon1, lat1)
132+
133+
# convert to other coords?:
134+
# https://gis.stackexchange.com/questions/78838/converting-projected-coordinates-to-lat-lon-using-python
135+
# https://openmaptiles.com/coordinate-systems/
136+
# https://ocefpaf.github.io/python4oceanographers/blog/2013/12/16/utm/
137+
# Web Mercator projection (EPSG:3857)
138+
# convert to wmp
139+
x0_wmp, y0_wmp = pyproj.transform(inProj, outProj, lon0, lat0)
140+
x1_wmp, y1_wmp = pyproj.transform(inProj, outProj, lon1, lat1)
141+
142+
# create data array
143+
out_arr = [lon0, lat0, lon1, lat1, x0_wmp, y0_wmp, x1_wmp, y1_wmp]
144+
145+
return out_arr, poly_geo
146+
27147

28148
###############################################################################
29-
def add_geo_coords_to_df(df_, create_geojson=False,
30-
inProj_str='epsg:4326', outProj_str='epsg:3857',
31-
verbose=False):
149+
def add_geo_coords_to_df(df_, inProj_str='epsg:4326', outProj_str='epsg:3857',
150+
create_geojson=False, verbose=False):
32151
'''Determine geo coords (latlon + wmp) of the bounding boxes'''
33152

34-
35-
# Placeholder
36-
return df_, []
37-
153+
"""
154+
Get geo coords for SIMRSWN dataframe containing pixel coords.
155+
156+
Arguments
157+
---------
158+
df_ : pandas dataframe
159+
Dataframe containing bounding box predictions
160+
inProj_str : str
161+
Projection of input image. Defaults to ``'epsg:4326'`` (WGS84)
162+
outProj_str : str
163+
Desired projection of bounding box coordinate.
164+
Defaults to ``'epsg:3857'`` (WGS84 Web Mercator for mapping).
165+
create_geojson : boolean
166+
Switch to output a geojson. If False, return [].
167+
Defaults to ``False``.
168+
verbose : boolean
169+
Switch to print values to screen. Defaults to ``False``.
170+
171+
Returns
172+
-------
173+
df_, json : tuple
174+
df_ is the updated dataframe with geo coords.
175+
json is the optional GeoJSON corresponding to the boxes.
176+
"""
177+
178+
t0 = time.time()
179+
print("Adding geo coords...")
180+
181+
# first, create a dictionary of geo transforms
182+
raster_locs = np.unique(df_['Image_Path'].values)
183+
geo_dict = {}
184+
for raster_loc in raster_locs:
185+
# raster geo transform
186+
gdal_geomTransform = gdal.Open(raster_loc).GetGeoTransform()
187+
affineObject = rio.open(raster_loc).transform
188+
if verbose:
189+
print("raster_loc:", raster_loc)
190+
print("gdal_geomTransform:", gdal_geomTransform)
191+
print("affineObject:", affineObject)
192+
geo_dict[raster_loc] = [gdal_geomTransform, affineObject]
193+
194+
# iterate through dataframe
195+
# columns = ['geometry']
196+
out_arr_json = []
197+
out_arr = []
198+
for idx, row in df_.iterrows():
199+
# get transform
200+
[gdal_geomTransform, affineObject] = geo_dict[row['Image_Path']]
201+
if verbose:
202+
print(idx, " row['Xmin_Glob'], row['Ymin_Glob'],",
203+
"row['Xmax_Glob'], row['Ymax_Glob']",
204+
row['Xmin_Glob'], row['Ymin_Glob'],
205+
row['Xmax_Glob'], row['Ymax_Glob'])
206+
207+
out_arr_row, poly_geo = get_row_geo_coords(
208+
row, affineObject=affineObject,
209+
gdal_geomTransform=gdal_geomTransform,
210+
inProj_str=inProj_str, outProj_str=outProj_str,
211+
verbose=verbose)
212+
out_arr.append(out_arr_row)
213+
if create_geojson:
214+
out_arr_json.append(poly_geo)
215+
216+
# update dataframe
217+
# [lon0, lat0, lon1, lat1, x0_wmp, y0_wmp, x1_wmp, y1_wmp]
218+
out_arr = np.array(out_arr)
219+
df_['lon0'] = out_arr[:, 0]
220+
df_['lat0'] = out_arr[:, 1]
221+
df_['lon1'] = out_arr[:, 2]
222+
df_['lat1'] = out_arr[:, 3]
223+
df_['x0_wmp'] = out_arr[:, 4]
224+
df_['y0_wmp'] = out_arr[:, 5]
225+
df_['x1_wmp'] = out_arr[:, 6]
226+
df_['y1_wmp'] = out_arr[:, 7]
227+
228+
# geodataframe if desired
229+
# https://gis.stackexchange.com/questions/174159/convert-a-pandas-dataframe-to-a-geodataframe
230+
if create_geojson:
231+
crs_init = {'init': inProj_str}
232+
df_json = pd.DataFrame(out_arr_json, columns=['geometry'])
233+
# add important columns to df_json
234+
df_json['category'] = df_['Category'].values
235+
df_json['prob'] = df_['Prob'].values
236+
gdf = gpd.GeoDataFrame(df_json, crs=crs_init, geometry=out_arr_json)
237+
json_out = gdf
238+
# json_out = gdf.to_json()
239+
else:
240+
json_out = []
241+
242+
t1 = time.time()
243+
print("Time to add geo coords to df:", t1 - t0, "seconds")
244+
return df_, json_out

0 commit comments

Comments
 (0)