Skip to content

Commit d34b59e

Browse files
committed
fixing geodesic buffer for new proj lib
1 parent 2d9d5fb commit d34b59e

File tree

2 files changed

+15
-14
lines changed

2 files changed

+15
-14
lines changed

ost/Project.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -366,7 +366,7 @@ def download(self, inventory_df, mirror=None, concurrent=2,
366366
download_df = inventory_df[inventory_df.download_path == 'None']
367367

368368
# to download or not ot download - that is here the question
369-
if not download_df.any().any():
369+
if download_df.empty:
370370
logger.info('All scenes are ready for being processed.')
371371
else:
372372
logger.info('One or more scene(s) need(s) to be downloaded.')

ost/helpers/vector.py

Lines changed: 14 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -171,7 +171,7 @@ def reproject_geometry(geom, inproj4, out_epsg):
171171
return geom
172172

173173

174-
def geodesic_point_buffer(lat, lon, meters, envelope=False):
174+
def geodesic_point_buffer(lon, lat, meters, envelope=False):
175175
"""
176176
177177
:param lat:
@@ -181,23 +181,24 @@ def geodesic_point_buffer(lat, lon, meters, envelope=False):
181181
:return:
182182
"""
183183

184-
# get WGS 84 proj
185-
proj_wgs84 = pyproj.Proj('epsg:4326')
186-
187-
# Azimuthal equidistant projection
188-
aeqd_proj = '+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0'
189-
project = partial(
190-
pyproj.transform,
191-
pyproj.Proj(aeqd_proj.format(lat=lat, lon=lon)),
192-
proj_wgs84
184+
proj_crs = ProjectedCRS(
185+
conversion=AzumuthalEquidistantConversion(float(lat), float(lon))
193186
)
194187

188+
proj_wgs84 = pyproj.Proj('EPSG:4326')
189+
190+
Trans = pyproj.Transformer.from_proj(
191+
proj_crs,
192+
proj_wgs84,
193+
always_xy=True
194+
).transform
195+
195196
buf = Point(0, 0).buffer(meters) # distance in metres
196197

197198
if envelope is True:
198-
geom = Polygon(transform(project, buf).exterior.coords[:]).envelope
199+
geom = Polygon(transform(Trans, buf).exterior.coords[:]).envelope
199200
else:
200-
geom = Polygon(transform(project, buf).exterior.coords[:])
201+
geom = Polygon(transform(Trans, buf).exterior.coords[:])
201202

202203
return geom.wkt
203204

@@ -240,7 +241,7 @@ def latlon_to_wkt(
240241
aoi_wkt = aoi_geom.to_wkt()
241242

242243
elif buffer_meter:
243-
aoi_wkt = geodesic_point_buffer(lat, lon, buffer_meter, envelope)
244+
aoi_wkt = geodesic_point_buffer(lon, lat, buffer_meter, envelope)
244245

245246
return aoi_wkt
246247

0 commit comments

Comments
 (0)