Efficient AI4EO OpenSource framework
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

341 lines
12KB

  1. """Slippy Map Tiles.
  2. See: https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
  3. """
  4. import io
  5. import os
  6. import re
  7. import glob
  8. import warnings
  9. import numpy as np
  10. from PIL import Image
  11. from rasterio import open as rasterio_open
  12. import cv2
  13. import json
  14. import psycopg2
  15. import rasterio
  16. import mercantile
  17. import supermercado
  18. warnings.simplefilter("ignore", UserWarning) # To prevent rasterio NotGeoreferencedWarning
  19. def tile_pixel_to_location(tile, dx, dy):
  20. """Converts a pixel in a tile to lon/lat coordinates."""
  21. assert 0 <= dx <= 1 and 0 <= dy <= 1, "x and y offsets must be in [0, 1]"
  22. w, s, e, n = mercantile.bounds(tile)
  23. def lerp(a, b, c):
  24. return a + c * (b - a)
  25. return lerp(w, e, dx), lerp(s, n, dy) # lon, lat
  26. def tiles_from_csv(path, xyz=True, extra_columns=False):
  27. """Retrieve tiles from a line-delimited csv file."""
  28. assert os.path.isfile(os.path.expanduser(path)), "'{}' seems not a valid CSV file".format(path)
  29. with open(os.path.expanduser(path)) as fp:
  30. for row in fp:
  31. row = row.replace("\n", "")
  32. if not row:
  33. continue
  34. row = re.split(",|\t", row) # use either comma or tab as separator
  35. if xyz:
  36. assert len(row) >= 3, "Invalid Cover"
  37. if not extra_columns or len(row) == 3:
  38. yield mercantile.Tile(int(row[0]), int(row[1]), int(row[2]))
  39. else:
  40. yield [mercantile.Tile(int(row[0]), int(row[1]), int(row[2])), *map(float, row[3:])]
  41. if not xyz:
  42. assert len(row) >= 1, "Invalid Cover"
  43. if not extra_columns:
  44. yield row[0]
  45. else:
  46. yield [row[0], *map(float, row[1:])]
  47. def tiles_from_dir(root, cover=None, xyz=True, xyz_path=False):
  48. """Loads files from an on-disk dir."""
  49. root = os.path.expanduser(root)
  50. if xyz is True:
  51. paths = glob.glob(os.path.join(root, "[0-9]*/[0-9]*/[0-9]*.*"))
  52. for path in paths:
  53. tile_xyz = re.match(os.path.join(root, "(?P<z>[0-9]+)/(?P<x>[0-9]+)/(?P<y>[0-9]+).+"), path)
  54. if not tile_xyz:
  55. continue
  56. tile = mercantile.Tile(int(tile_xyz["x"]), int(tile_xyz["y"]), int(tile_xyz["z"]))
  57. if cover is not None and tile not in cover:
  58. continue
  59. if xyz_path is True:
  60. yield tile, path
  61. else:
  62. yield tile
  63. else:
  64. paths = glob.glob(root, "**/*.*", recursive=True)
  65. for path in paths:
  66. return path
  67. def tile_from_xyz(root, x, y, z):
  68. """Retrieve a single tile from a slippy map dir."""
  69. path = glob.glob(os.path.join(os.path.expanduser(root), str(z), str(x), str(y) + ".*"))
  70. if not path:
  71. return None
  72. assert len(path) == 1, "ambiguous tile path"
  73. return mercantile.Tile(x, y, z), path[0]
  74. def tile_bbox(tile, mercator=False):
  75. if isinstance(tile, mercantile.Tile):
  76. if mercator:
  77. return mercantile.xy_bounds(tile) # EPSG:3857
  78. else:
  79. return mercantile.bounds(tile) # EPSG:4326
  80. else:
  81. with open(rasterio_open(tile)) as r:
  82. if mercator:
  83. w, s, e, n = r.bounds
  84. w, s = mercantile.xy(w, s)
  85. e, n = mercantile.xy(e, n)
  86. return w, s, e, n # EPSG:3857
  87. else:
  88. return r.bounds # EPSG:4326
  89. assert False, "Unable to open tile"
  90. def tiles_to_geojson(tiles, union=True):
  91. """Convert tiles to their footprint GeoJSON."""
  92. first = True
  93. geojson = '{"type":"FeatureCollection","features":['
  94. if union: # smaller tiles union geometries (but losing properties)
  95. tiles = [str(tile.z) + "-" + str(tile.x) + "-" + str(tile.y) + "\n" for tile in tiles]
  96. for feature in supermercado.uniontiles.union(tiles, True):
  97. geojson += json.dumps(feature) if first else "," + json.dumps(feature)
  98. first = False
  99. else: # keep each tile geometry and properties (but fat)
  100. for tile in tiles:
  101. prop = '"properties":{{"x":{},"y":{},"z":{}}}'.format(tile.x, tile.y, tile.z)
  102. geom = '"geometry":{}'.format(json.dumps(mercantile.feature(tile, precision=6)["geometry"]))
  103. geojson += '{}{{"type":"Feature",{},{}}}'.format("," if not first else "", geom, prop)
  104. first = False
  105. geojson += "]}"
  106. return geojson
  107. def tiles_to_granules(tiles, pg):
  108. """Retrieve Intersecting Sentinel Granules from tiles."""
  109. conn = psycopg2.connect(pg)
  110. db = conn.cursor()
  111. assert db
  112. granules = set()
  113. tiles = [str(tile.z) + "-" + str(tile.x) + "-" + str(tile.y) + "\n" for tile in tiles]
  114. for feature in supermercado.uniontiles.union(tiles, True):
  115. geom = json.dumps(feature["geometry"])
  116. query = """SELECT id FROM neo.s2_granules
  117. WHERE ST_Intersects(geom, ST_SetSRID(ST_GeomFromGeoJSON('{}'), 4326))""".format(
  118. geom
  119. )
  120. db.execute(query)
  121. granules.update([str(granule[0]) for granule in db.fetchall()])
  122. return granules
  123. def tile_image_from_file(path, bands=None, force_rgb=False):
  124. """Return a multiband image numpy array, from an image file path, or None."""
  125. try:
  126. if path[-3:] == "png" and force_rgb: # PIL PNG Color Palette handling
  127. return np.array(Image.open(os.path.expanduser(path)).convert("RGB"))
  128. elif path[-3:] == "png":
  129. return np.array(Image.open(os.path.expanduser(path)))
  130. else:
  131. raster = rasterio_open(os.path.expanduser(path))
  132. except:
  133. return None
  134. image = None
  135. for i in raster.indexes if bands is None else bands:
  136. data_band = raster.read(i)
  137. data_band = data_band.reshape(data_band.shape[0], data_band.shape[1], 1) # H,W -> H,W,C
  138. image = np.concatenate((image, data_band), axis=2) if image is not None else data_band
  139. assert image is not None, "Unable to open {}".format(path)
  140. return image
  141. def tile_image_to_file(root, tile, image):
  142. """ Write an image tile on disk. """
  143. H, W, C = image.shape
  144. root = os.path.expanduser(root)
  145. path = os.path.join(root, str(tile.z), str(tile.x)) if isinstance(tile, mercantile.Tile) else root
  146. os.makedirs(path, exist_ok=True)
  147. if C == 1:
  148. ext = "png"
  149. elif C == 3:
  150. ext = "jpg"
  151. else:
  152. ext = "tif"
  153. if isinstance(tile, mercantile.Tile):
  154. path = os.path.join(path, "{}.{}".format(str(tile.y), ext))
  155. else:
  156. path = os.path.join(path, "{}.{}".format(tile, ext))
  157. try:
  158. if C == 1:
  159. Image.fromarray(image.reshape(H, W), mode="L").save(path)
  160. elif C == 3:
  161. cv2.imwrite(path, cv2.cvtColor(image, cv2.COLOR_RGB2BGR))
  162. else:
  163. rasterio.open(path, "w", driver="GTiff", height=H, width=W, count=C, dtype=image.dtype).write(
  164. np.moveaxis(image, 2, 0) # H,W,C -> C,H,W
  165. )
  166. except:
  167. assert False, "Unable to write {}".format(path)
  168. def tile_label_from_file(path, silent=True):
  169. """Return a numpy array, from a label file path, or None."""
  170. try:
  171. return np.array(Image.open(os.path.expanduser(path))).astype(int)
  172. except:
  173. assert silent, "Unable to open existing label: {}".format(path)
  174. def tile_label_to_file(root, tile, palette, transparency, label, append=False, margin=0):
  175. """ Write a label (or a mask) tile on disk. """
  176. root = os.path.expanduser(root)
  177. dir_path = os.path.join(root, str(tile.z), str(tile.x)) if isinstance(tile, mercantile.Tile) else root
  178. path = os.path.join(dir_path, "{}.png".format(str(tile.y)))
  179. if len(label.shape) == 3: # H,W,C -> H,W
  180. assert label.shape[2] == 1
  181. label = label.reshape((label.shape[0], label.shape[1]))
  182. if append and os.path.isfile(path):
  183. previous = tile_label_from_file(path, silent=False)
  184. label = np.uint8(np.maximum(previous, label))
  185. else:
  186. os.makedirs(dir_path, exist_ok=True)
  187. try:
  188. out = Image.fromarray(label, mode="P")
  189. out.putpalette(palette)
  190. if transparency is not None:
  191. out.save(path, optimize=True, transparency=transparency)
  192. else:
  193. out.save(path, optimize=True)
  194. except:
  195. assert False, "Unable to write {}".format(path)
  196. def tile_image_from_url(requests_session, url, timeout=10):
  197. """Fetch a tile image using HTTP, and return it or None """
  198. try:
  199. resp = requests_session.get(url, timeout=timeout)
  200. resp.raise_for_status()
  201. image = np.fromstring(io.BytesIO(resp.content).read(), np.uint8)
  202. return cv2.cvtColor(cv2.imdecode(image, cv2.IMREAD_ANYCOLOR), cv2.COLOR_BGR2RGB)
  203. except Exception:
  204. return None
  205. def tile_is_neighboured(tile, tiles):
  206. """Check if a tile is surrounded by others tiles"""
  207. tiles = dict(tiles)
  208. try:
  209. # 3x3 matrix (upper, center, bottom) x (left, center, right)
  210. tiles[mercantile.Tile(x=int(tile.x) - 1, y=int(tile.y) - 1, z=int(tile.z))] # ul
  211. tiles[mercantile.Tile(x=int(tile.x) + 0, y=int(tile.y) - 1, z=int(tile.z))] # uc
  212. tiles[mercantile.Tile(x=int(tile.x) + 1, y=int(tile.y) - 1, z=int(tile.z))] # ur
  213. tiles[mercantile.Tile(x=int(tile.x) - 1, y=int(tile.y) + 0, z=int(tile.z))] # cl
  214. tiles[mercantile.Tile(x=int(tile.x) + 1, y=int(tile.y) + 0, z=int(tile.z))] # cr
  215. tiles[mercantile.Tile(x=int(tile.x) - 1, y=int(tile.y) + 1, z=int(tile.z))] # bl
  216. tiles[mercantile.Tile(x=int(tile.x) + 0, y=int(tile.y) + 1, z=int(tile.z))] # bc
  217. tiles[mercantile.Tile(x=int(tile.x) + 1, y=int(tile.y) + 1, z=int(tile.z))] # br
  218. except KeyError:
  219. return False
  220. return True
  221. def tile_image_buffer(tile, tiles, bands):
  222. """Buffers a tile image adding borders on all sides based on adjacent tile, or zeros padded if not possible."""
  223. def tile_image_neighbour(tile, dx, dy, tiles, bands):
  224. """Retrieves neighbour tile image if exists."""
  225. try:
  226. path = tiles[mercantile.Tile(x=int(tile.x) + dx, y=int(tile.y) + dy, z=int(tile.z))]
  227. except KeyError:
  228. return None
  229. return tile_image_from_file(path, bands)
  230. tiles = dict(tiles)
  231. # 3x3 matrix (upper, center, bottom) x (left, center, right)
  232. ul = tile_image_neighbour(tile, -1, -1, tiles, bands)
  233. uc = tile_image_neighbour(tile, +0, -1, tiles, bands)
  234. ur = tile_image_neighbour(tile, +1, -1, tiles, bands)
  235. cl = tile_image_neighbour(tile, -1, +0, tiles, bands)
  236. cc = tile_image_neighbour(tile, +0, +0, tiles, bands)
  237. cr = tile_image_neighbour(tile, +1, +0, tiles, bands)
  238. bl = tile_image_neighbour(tile, -1, +1, tiles, bands)
  239. bc = tile_image_neighbour(tile, +0, +1, tiles, bands)
  240. br = tile_image_neighbour(tile, +1, +1, tiles, bands)
  241. b = len(bands)
  242. ts = cc.shape[1]
  243. o = int(ts / 4)
  244. oo = o * 2
  245. img = np.zeros((ts + oo, ts + oo, len(bands))).astype(np.uint8)
  246. # fmt:off
  247. img[0:o, 0:o, :] = ul[-o:ts, -o:ts, :] if ul is not None else np.zeros((o, o, b)).astype(np.uint8)
  248. img[0:o, o:ts+o, :] = uc[-o:ts, 0:ts, :] if uc is not None else np.zeros((o, ts, b)).astype(np.uint8)
  249. img[0:o, ts+o:ts+oo, :] = ur[-o:ts, 0:o, :] if ur is not None else np.zeros((o, o, b)).astype(np.uint8)
  250. img[o:ts+o, 0:o, :] = cl[0:ts, -o:ts, :] if cl is not None else np.zeros((ts, o, b)).astype(np.uint8)
  251. img[o:ts+o, o:ts+o, :] = cc if cc is not None else np.zeros((ts, ts, b)).astype(np.uint8)
  252. img[o:ts+o, ts+o:ts+oo, :] = cr[0:ts, 0:o, :] if cr is not None else np.zeros((ts, o, b)).astype(np.uint8)
  253. img[ts+o:ts+oo, 0:o, :] = bl[0:o, -o:ts, :] if bl is not None else np.zeros((o, o, b)).astype(np.uint8)
  254. img[ts+o:ts+oo, o:ts+o, :] = bc[0:o, 0:ts, :] if bc is not None else np.zeros((o, ts, b)).astype(np.uint8)
  255. img[ts+o:ts+oo, ts+o:ts+oo, :] = br[0:o, 0:o, :] if br is not None else np.zeros((o, o, b)).astype(np.uint8)
  256. # fmt:on
  257. return img