Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove rasterio dep on pyproj #1182

Merged
merged 4 commits into from
May 30, 2023
Merged

Remove rasterio dep on pyproj #1182

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b739cd9
Remove rasterio dep on pyproj
banesullivan May 29, 2023
5763b4b
Great circle estimation when pyproj unavailable
banesullivan May 30, 2023
43d3efe
Remove unused imports
banesullivan May 30, 2023
ee6a6c0
Cleanup
banesullivan May 30, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 28 additions & 13 deletions large_image/tilesource/geo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,14 +1,18 @@
from urllib.parse import urlencode, urlparse

import pyproj # TODO: import issues

from large_image.cache_util import CacheProperties, methodcache
from large_image.constants import SourcePriority, TileInputUnits
from large_image.exceptions import TileSourceError

from .base import FileTileSource
from .utilities import JSONDict, getPaletteColors

try:
import pyproj
has_pyproj = True
except Exception:
has_pyproj = False

# Inform the tile source cache about the potential size of this tile source
CacheProperties['tilesource']['itemExpectedSize'] = max(
CacheProperties['tilesource']['itemExpectedSize'],
Expand All @@ -18,8 +22,9 @@
ProjUnitsAcrossLevel0 = {}
ProjUnitsAcrossLevel0_MaxSize = 100

InitPrefix = '+init='
NeededInitPrefix = '' if int(pyproj.proj_version_str.split('.')[0]) >= 6 else InitPrefix
InitPrefix = ''
if has_pyproj:
NeededInitPrefix = '+init=' if int(pyproj.proj_version_str.split('.')[0]) < 6 else InitPrefix


def make_vsi(url: str, **options):
Expand Down Expand Up @@ -225,15 +230,25 @@ def getPixelSizeInMeters(self):
bounds = self.getBounds(NeededInitPrefix + 'epsg:4326')
if not bounds:
return
geod = pyproj.Geod(ellps='WGS84')
az12, az21, s1 = geod.inv(bounds['ul']['x'], bounds['ul']['y'],
bounds['ur']['x'], bounds['ur']['y'])
az12, az21, s2 = geod.inv(bounds['ur']['x'], bounds['ur']['y'],
bounds['lr']['x'], bounds['lr']['y'])
az12, az21, s3 = geod.inv(bounds['lr']['x'], bounds['lr']['y'],
bounds['ll']['x'], bounds['ll']['y'])
az12, az21, s4 = geod.inv(bounds['ll']['x'], bounds['ll']['y'],
bounds['ul']['x'], bounds['ul']['y'])
if has_pyproj:
geod = pyproj.Geod(ellps='WGS84')
computer = geod.inv
else:
# Estimate based on great-cirlce distance
def computer(lon1, lat1, lon2, lat2):
from math import acos, cos, radians, sin
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
return None, None, 6.378e+6 * (
acos(sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon1 - lon2))
)
_, _, s1 = computer(bounds['ul']['x'], bounds['ul']['y'],
bounds['ur']['x'], bounds['ur']['y'])
_, _, s2 = computer(bounds['ur']['x'], bounds['ur']['y'],
bounds['lr']['x'], bounds['lr']['y'])
_, _, s3 = computer(bounds['lr']['x'], bounds['lr']['y'],
bounds['ll']['x'], bounds['ll']['y'])
_, _, s4 = computer(bounds['ll']['x'], bounds['ll']['y'],
bounds['ul']['x'], bounds['ul']['y'])
return (s1 + s2 + s3 + s4) / (self.sourceSizeX * 2 + self.sourceSizeY * 2)

def getNativeMagnification(self):
Expand Down
59 changes: 32 additions & 27 deletions sources/rasterio/large_image_source_rasterio/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,7 @@
import PIL.Image
import rasterio as rio
from affine import Affine
from pyproj import CRS, Transformer
from pyproj.exceptions import CRSError
from rasterio import warp
from rasterio.enums import ColorInterp, Resampling
from rasterio.errors import RasterioIOError

Expand Down Expand Up @@ -58,6 +57,16 @@
warnings.filterwarnings('ignore', category=rio.errors.NotGeoreferencedWarning, module='rasterio')


def make_crs(projection):
if isinstance(projection, str):
return rio.CRS.from_string(projection)
if isinstance(projection, dict):
return rio.CRS.from_dict(projection)
if isinstance(projection, int):
return rio.CRS.from_string(f'EPSG:{projection}')
return rio.CRS(projection)


class RasterioFileTileSource(GDALBaseFileTileSource, metaclass=LruCacheMetaclass):
"""Provides tile access to geospatial files."""

Expand All @@ -70,7 +79,7 @@ def __init__(self, path, projection=None, unitsPerPixel=None, **kwargs):
See the base class for other available parameters.

:param path: a filesystem path for the tile source.
:param projection: None to use pixel space, otherwise a crs compatible with pyproj.
:param projection: None to use pixel space, otherwise a crs compatible with rasterio's CRS.
:param unitsPerPixel: The size of a pixel at the 0 tile size.
Ignored if the projection is None. For projections, None uses the default,
which is the distance between (-180,0) and (180,0) in EPSG:4326 converted to the
Expand Down Expand Up @@ -111,7 +120,7 @@ def __init__(self, path, projection=None, unitsPerPixel=None, **kwargs):
self._bounds = {}
self.tileWidth = self.tileSize
self.tileHeight = self.tileSize
self.projection = CRS(projection) if projection else None
self.projection = make_crs(projection) if projection else None

# get width and height parameters
with self._getDatasetLock:
Expand Down Expand Up @@ -211,7 +220,7 @@ def _initWithProjection(self, unitsPerPixel=None):

:param unitsPerPixel: optional default to None
"""
srcCrs = CRS(4326)
srcCrs = make_crs(4326)
# Since we already converted to bytes decoding is safe here
dstCrs = self.projection
if dstCrs.is_geographic:
Expand All @@ -230,8 +239,8 @@ def _initWithProjection(self, unitsPerPixel=None):
# If unitsPerPixel is not specified, the horizontal distance
# between -180,0 and +180,0 is used. Some projections (such as
# stereographic) will fail in this case; they must have a unitsPerPixel specified.
equator = Transformer.from_crs(srcCrs, dstCrs, always_xy=True)
east, west = equator.itransform([(-180, 0), (180, 0)])
east, _ = warp.transform(srcCrs, dstCrs, [-180,], [0,])
west, _ = warp.transform(srcCrs, dstCrs, [180,], [0,])
self.unitsAcrossLevel0 = abs(east[0] - west[0])
if not self.unitsAcrossLevel0:
raise TileSourceError(
Expand Down Expand Up @@ -293,7 +302,7 @@ def getCrs(self):
hasTransform = self.dataset.transform != Affine.identity()
isNitf = self.dataset.driver.lower() in {'NITF'}
if not crs and (hasTransform or isNitf):
crs = CRS(4326)
crs = make_crs(4326)

return crs

Expand Down Expand Up @@ -326,7 +335,7 @@ def getBounds(self, crs=None, **kwargs):
crs = kwargs.get('srs')

# read the crs as a crs if needed
dstCrs = CRS(crs) if crs else None
dstCrs = make_crs(crs) if crs else None
strDstCrs = 'none' if dstCrs is None else dstCrs.to_string()

# exit if it's already set
Expand Down Expand Up @@ -370,8 +379,7 @@ def getBounds(self, crs=None, **kwargs):
# set the vertical bounds
# some projection system don't cover the poles so we need to adapt
# the values of ybounds accordingly
transformer = Transformer.from_crs(4326, dstCrs, always_xy=True)
has_poles = transformer.transform(0, 90)[1] != float('inf')
has_poles = warp.transform(4326, dstCrs, [0], [90])[1][0] != float('inf')
yBounds = 90 if has_poles else 89.999999

# for each corner fix the latitude within -yBounds yBounds
Expand All @@ -394,9 +402,8 @@ def getBounds(self, crs=None, **kwargs):
# reproject the pts in the destination coordinate system if necessary
needProjection = dstCrs and dstCrs != srcCrs
if needProjection:
transformer = Transformer.from_crs(srcCrs, dstCrs, always_xy=True)
for pt in bounds.values():
pt['x'], pt['y'] = transformer.transform(pt['x'], pt['y'])
[pt['x']], [pt['y']] = warp.transform(srcCrs, dstCrs, [pt['x']], [pt['y']])

# extract min max coordinates from the corners
ll = bounds['ll']['x'], bounds['ll']['y']
Expand Down Expand Up @@ -671,15 +678,14 @@ def _convertProjectionUnits(
if units.startswith('proj4:'):
# HACK to avoid `proj4:` prefixes with `WGS84`, etc.
units = units.split(':', 1)[1]
srcCrs = CRS(units)
srcCrs = make_crs(units)
dstCrs = self.projection # instance projection -- do not use the CRS native to the file
transformer = Transformer.from_crs(srcCrs, dstCrs, always_xy=True)
pleft, ptop = transformer.transform(
right if left is None else left,
bottom if top is None else top)
pright, pbottom = transformer.transform(
left if right is None else right,
top if bottom is None else bottom)
[pleft], [ptop] = warp.transform(srcCrs, dstCrs,
[right if left is None else left],
[bottom if top is None else top])
[pright], [pbottom] = warp.transform(srcCrs, dstCrs,
[left if right is None else right],
[top if bottom is None else bottom])
units = 'projection'

# set the corner value in pixel coordinates if the coordinate was initially
Expand Down Expand Up @@ -730,8 +736,8 @@ def _getRegionBounds(

# check if the units is a string or projection material
isProj = False
with suppress(CRSError):
isProj = CRS(units) is not None
with suppress(rio.errors.CRSError):
isProj = make_crs(units) is not None

# convert the coordinates if a projection exist
if isUnits and isProj:
Expand Down Expand Up @@ -821,12 +827,11 @@ def toNativePixelCoordinates(self, x, y, crs=None, roundResults=True):

:return: (x, y) the pixel coordinate.
"""
srcCrs = self.projection if crs is None else CRS(crs)
srcCrs = self.projection if crs is None else make_crs(crs)

# convert to the native projection
dstCrs = CRS(self.getCrs())
transformer = Transformer.from_crs(srcCrs, dstCrs, always_xy=True)
px, py = transformer.transform(x, y)
dstCrs = make_crs(self.getCrs())
[px], [py] = warp.transform(srcCrs, dstCrs, [x], [y])

# convert to native pixel coordinates
af = self._getAffine()
Expand Down