diff --git a/large_image/tilesource/geo.py b/large_image/tilesource/geo.py index 27e8d2b27..289a0a6ba 100644 --- a/large_image/tilesource/geo.py +++ b/large_image/tilesource/geo.py @@ -1,7 +1,5 @@ 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 @@ -9,6 +7,12 @@ 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'], @@ -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): @@ -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): diff --git a/sources/rasterio/large_image_source_rasterio/__init__.py b/sources/rasterio/large_image_source_rasterio/__init__.py index 5eba97fee..533b6dec5 100644 --- a/sources/rasterio/large_image_source_rasterio/__init__.py +++ b/sources/rasterio/large_image_source_rasterio/__init__.py @@ -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 @@ -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.""" @@ -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 @@ -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: @@ -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: @@ -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( @@ -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 @@ -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 @@ -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 @@ -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'] @@ -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 @@ -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: @@ -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()