/
utils.py
191 lines (148 loc) · 5.19 KB
/
utils.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
"""rio_cogeo.utils: Utility functions."""
import math
import warnings
from typing import Dict, Tuple
import mercantile
from rasterio.crs import CRS
from rasterio.enums import ColorInterp, MaskFlags
from rasterio.enums import Resampling as ResamplingEnums
from rasterio.transform import Affine
from rasterio.warp import calculate_default_transform, transform_bounds
from supermercado.burntiles import tile_extrema
from rio_cogeo.errors import DeprecationWarning
def _meters_per_pixel(zoom, lat=0.0, tilesize=256):
"""
Return the pixel resolution for a given mercator tile zoom and lattitude.
Parameters
----------
zoom: int
Mercator zoom level
lat: float, optional
Latitude in decimal degree (default: 0)
tilesize: int, optional
Mercator tile size (default: 256).
Returns
-------
Pixel resolution in meters
"""
return (math.cos(lat * math.pi / 180.0) * 2 * math.pi * 6378137) / (
tilesize * 2 ** zoom
)
def zoom_for_pixelsize(pixel_size, max_z=24, tilesize=256):
"""
Get mercator zoom level corresponding to a pixel resolution.
Freely adapted from
https://github.com/OSGeo/gdal/blob/b0dfc591929ebdbccd8a0557510c5efdb893b852/gdal/swig/python/scripts/gdal2tiles.py#L294
Parameters
----------
pixel_size: float
Pixel size
max_z: int, optional (default: 24)
Max mercator zoom level allowed
tilesize: int, optional
Mercator tile size (default: 256).
Returns
-------
Mercator zoom level corresponding to the pixel resolution
"""
for z in range(max_z):
if pixel_size > _meters_per_pixel(z, 0, tilesize=tilesize):
return max(0, z - 1) # We don't want to scale up
return max_z - 1
def get_zooms(src_dst, lat=0.0, tilesize=256) -> Tuple[int, int]:
"""
Calculate raster max zoom level.
Parameters
----------
src: rasterio.io.DatasetReader
Rasterio io.DatasetReader object
lat: float, optional
Center latitude of the dataset. This is only needed in case you want to
apply latitude correction factor to ensure consitent maximum zoom level
(default: 0.0).
tilesize: int, optional
Mercator tile size (default: 256).
Returns
-------
max_zoom: int
Max zoom level.
"""
dst_affine, w, h = calculate_default_transform(
src_dst.crs, "epsg:3857", src_dst.width, src_dst.height, *src_dst.bounds
)
native_resolution = max(abs(dst_affine[0]), abs(dst_affine[4]))
# Correction factor for web-mercator projection latitude distortion
latitude_correction_factor = math.cos(math.radians(lat))
corrected_resolution = native_resolution * latitude_correction_factor
max_zoom = zoom_for_pixelsize(corrected_resolution, tilesize=tilesize)
ovr_resolution = corrected_resolution * max(h, w) / tilesize
min_zoom = zoom_for_pixelsize(ovr_resolution, tilesize=tilesize)
return (min_zoom, max_zoom)
def get_maximum_overview_level(src_dst, minsize=512):
"""
Calculate the maximum overview level.
Attributes
----------
src_dst : rasterio.io.DatasetReader
Rasterio io.DatasetReader object.
minsize : int (default: 512)
Minimum overview size.
Returns
-------
nlevel: int
overview level.
"""
warnings.warn(
"rio_cogeo.utils.get_maximum_overview_level will be removed in version 2.0",
DeprecationWarning,
)
nlevel = 0
overview = 1
while min(src_dst.width // overview, src_dst.height // overview) > minsize:
overview *= 2
nlevel += 1
return nlevel
def has_alpha_band(src_dst):
"""Check for alpha band or mask in source."""
if (
any([MaskFlags.alpha in flags for flags in src_dst.mask_flag_enums])
or ColorInterp.alpha in src_dst.colorinterp
):
return True
return False
def has_mask_band(src_dst):
"""Check for mask band in source."""
if any([MaskFlags.per_dataset in flags for flags in src_dst.mask_flag_enums]):
return True
return False
def get_web_optimized_params(
src_dst,
tilesize=256,
latitude_adjustment: bool = True,
warp_resampling: str = "nearest",
grid_crs=CRS.from_epsg(3857),
) -> Dict:
"""Return VRT parameters for a WebOptimized COG."""
bounds = list(
transform_bounds(
src_dst.crs, CRS.from_epsg(4326), *src_dst.bounds, densify_pts=21
)
)
center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2]
lat = 0 if latitude_adjustment else center[1]
_, max_zoom = get_zooms(src_dst, lat=lat, tilesize=tilesize)
extrema = tile_extrema(bounds, max_zoom)
left, _, _, top = mercantile.xy_bounds(
extrema["x"]["min"], extrema["y"]["min"], max_zoom
)
vrt_res = _meters_per_pixel(max_zoom, 0, tilesize=tilesize)
vrt_transform = Affine(vrt_res, 0, left, 0, -vrt_res, top)
vrt_width = (extrema["x"]["max"] - extrema["x"]["min"]) * tilesize
vrt_height = (extrema["y"]["max"] - extrema["y"]["min"]) * tilesize
return dict(
crs=grid_crs,
transform=vrt_transform,
width=vrt_width,
height=vrt_height,
resampling=ResamplingEnums[warp_resampling],
)