· 5 years ago · Jun 11, 2020, 06:20 AM
1import concurrent.futures
2from collections import defaultdict
3from typing import Tuple, List, DefaultDict
4
5import imageio
6import mercantile
7import numpy as np
8import rasterio as rio
9import requests
10from affine import Affine
11from rasterio import transform, MemoryFile
12from rasterio.coords import BoundingBox
13from rasterio.errors import RasterioIOError
14from rasterio.windows import Window
15
16import boto3
17s3 = boto3.client("s3")
18
19
20def get_tiles_from_field_metadata(rgb_key: str, channel: str, zooms: int) -> List[List[int]]:
21 # Tiles per defined zoom level with mercantile API.
22 response = requests.get(f"https://jgau9m9qpf.execute-api.us-east-1.amazonaws.com/bondville"
23 f"/metadata/{rgb_key}/{channel}/")
24 west, south, east, north = response.json()['bounds']
25 tiles = list(mercantile.tiles(west, south, east, north, zooms))
26 final_tiles = [[tile.x, tile.y, tile.z] for tile in tiles]
27 return final_tiles
28
29
30def get_resulting_image_shape(rgb_key: str, channel: str, zooms: int,
31 tile_size: tuple = (256, 256)) -> Tuple[int, int]:
32 # Final resulting image shape (height, width).
33 tiles = get_tiles_from_field_metadata(rgb_key, channel, zooms)
34 combined_arrays_by_x = combine_image_data_by_x(tiles)
35
36 image_width, image_height = len(combined_arrays_by_x.keys()) * tile_size[0], None
37 for k in combined_arrays_by_x:
38 image_height = len(combined_arrays_by_x[k]) * tile_size[1]
39 break
40 return image_width, image_height
41
42
43def combine_image_data_by_x(tiles: List[List[int]]) -> DefaultDict[int, List[List[int]]]:
44 # Get dict of image array combined by x value.
45 initial_x = tiles[0][0]
46 combined_arrays_by_x = defaultdict(list)
47 for tile in tiles:
48 local_x = tile[0]
49 if local_x == initial_x:
50 xyz_array = [tile[0], tile[1], tile[2]]
51 combined_arrays_by_x[local_x].append(xyz_array)
52 else:
53 initial_x = local_x
54 xyz_array = [tile[0], tile[1], tile[2]]
55 combined_arrays_by_x[initial_x].append(xyz_array)
56 return combined_arrays_by_x
57
58
59def get_tile_bounds(tile: list) -> BoundingBox:
60 # Returns BoundingBox object for each given tile.
61 return BoundingBox(*mercantile.xy_bounds(tile))
62
63
64def transform_from_bounds(bounds, tile_size) -> Affine:
65 return transform.from_bounds(*bounds, *tile_size)
66
67
68def get_resulting_image_bounds(rgb_key: str, channel: str, zooms: int) -> BoundingBox:
69 # Final resulting image boundary (x_min, y_min, x_max, y_max).
70 tiles = get_tiles_from_field_metadata(rgb_key, channel, zooms)
71 x_min, y_min, x_max, y_max = [], [], [], []
72 for index, tile in enumerate(tiles):
73 bbox = get_tile_bounds(tiles[index])
74 x_min.append(bbox.left)
75 y_min.append(bbox.bottom)
76 x_max.append(bbox.right)
77 y_max.append(bbox.top)
78
79 return BoundingBox(min(x_min), min(y_min), max(x_max), max(y_max))
80
81
82def construct_url(rgb_key: str, x: int, y: int, z: int):
83 # Makes api calls to rgb tile endpoint.
84 return f"https://jgau9m9qpf.execute-api.us-east-1.amazonaws.com/bondville/" \
85 f"{rgb_key}/{x}/{y}/{z}.png?&r=red&g=green&b=blue" \
86 f"&r_range=%5B2048.0%2C8735.0%5D&g_range=%5B2048.0%2C8735.0%5D&b_range=%5B2048.0%2C8735.0%5D"
87
88
89def get_image_data_array(rgb_key: str, x: int, y: int, z: int):
90 # Construct tile data each time calling api with respective x, y, z coordinates.
91 image_data = imageio.imread(construct_url(rgb_key, x, y, z))
92 return image_data
93
94
95def write_image_for_one_tile(destination, image_data, image_width, image_height, g_transform,
96 no_data, window, count=3, mode='r+', driver='GTiff', crs='epsg:3857'):
97 # Write image data per tile basis.
98 with rio.open(destination, mode,
99 driver=driver,
100 width=image_width,
101 height=image_height,
102 crs=crs,
103 transform=g_transform,
104 nodata=no_data,
105 count=count,
106 dtype=image_data.dtype
107 ) as dst:
108 dst.write(np.transpose(np.moveaxis(image_data, [0, 2, 1], [1, 2, 0])), window=window)
109 # dst_profile = {
110 # 'driver': driver,
111 # 'count': count,
112 # 'height': image_height,
113 # 'width': image_width,
114 # 'dtype': image_data.dtype,
115 # 'nodata': no_data,
116 # 'crs': crs,
117 # 'transform': g_transform,
118 # }
119 # with MemoryFile() as memory_file:
120 # with memory_file.open(**dst_profile) as dst:
121 # dst.write(np.transpose(np.moveaxis(image_data, [0, 2, 1], [1, 2, 0])), window=window)
122
123
124def thread_function(d_file, row_offset: int, column_offset: int, tile: list, width: int,
125 height: int, g_transform: Affine, no_data_v: int):
126 # Thread function in order to pass to ThreadPullExecutor
127 image_data_ = get_image_data_array("rgb/flight/9BINMT6XB", tile[2], tile[0], tile[1])
128 tile_window = Window(column_offset, row_offset, 256, 256)
129 # write_image_for_one_tile(image_data_, width,
130 # height, g_transform, no_data_v,
131 # tile_window)
132 try:
133 write_image_for_one_tile(d_file, image_data_, width,
134 height, g_transform, no_data_v,
135 tile_window)
136 except RasterioIOError:
137 write_image_for_one_tile(d_file, image_data_, width,
138 height, g_transform, no_data_v,
139 tile_window, mode='w')
140
141
142if __name__ == "__main__":
143 resulting_tiles = get_tiles_from_field_metadata('flight/9BINMT6XB', 'red', 18)
144 res_image_width, res_image_height = get_resulting_image_shape('flight/9BINMT6XB', 'red', 18)
145 res_bounds = get_resulting_image_bounds('flight/9BINMT6XB', 'red', 18)
146 geo_transform = transform_from_bounds(res_bounds, (res_image_width, res_image_height))
147 tiles_map = combine_image_data_by_x(resulting_tiles)
148
149 # Please refer here to come with solution to "error setting certificate verify locations".
150 # https://github.com/mapbox/rasterio/issues/1289
151 # s3_destination = 'intelinair-infra'
152
153 s3_destination = 'temp.tif'
154 no_data_value = 0
155
156 col_off = -256
157
158 for key, x_tiles in tiles_map.items():
159 col_off += 256
160 row_off = -256
161
162 with concurrent.futures.ThreadPoolExecutor(max_workers=5) as executor:
163 futures = []
164 for x_tile in x_tiles:
165 row_off += 256
166 futures.append(executor.submit(thread_function, s3_destination, row_off, col_off,
167 x_tile, res_image_width, res_image_height,
168 geo_transform, no_data_value))
169 concurrent.futures.wait(futures)
170
171 # s3.upload_fileobj(mem_file, s3_destination, 'temp.tif')
172
173 # The below commented part is the sequential implementation of raster creation.
174 # for x_tile in x_tiles:
175 # row_off += 256
176 # image_data_ = get_image_data_array("rgb/flight/9BINMT6XB", x_tile[2],
177 # x_tile[0], x_tile[1])
178 # tile_window = Window(col_off, row_off, 256, 256)
179 # write_image_for_one_tile(s3_destination, image_data_, res_image_width,
180 # res_image_height, geo_transform, no_data_value,
181 # tile_window)
182 # try:
183 # write_image_for_one_tile(s3_destination, image_data_, res_image_width,
184 # res_image_height, geo_transform, no_data_value,
185 # tile_window)
186 # except RasterioIOError:
187 # write_image_for_one_tile(s3_destination, image_data_, res_image_width,
188 # res_image_height, geo_transform, no_data_value,
189 # tile_window, mode='w')
190
191# The below commented part is the raster generation other solution with numpy array vertical and
192# horizontal concatenation but with no geographic data applied.
193# Also the deletion of black pixels for the resulting image with addition of alpha channel.
194
195
196# def image_from_array(array: np.ndarray, file_name: str, extension: str, mode: str = 'RGB',
197# compress_level: int = 1):
198# # Convert array to PIL Image.
199#
200# # Make black pixels transparent, adding alpha channel to image, and converting all the pixels
201# # which rbg values are black to transparent.
202#
203# rgb_img = Image.fromarray(array, mode)
204# rgb_alpha_image = rgb_img.convert('RGBA')
205# pixel_data = list(rgb_alpha_image.getdata())
206# for i, pixel in enumerate(pixel_data):
207# if pixel[:3] == (0, 0, 0):
208# pixel_data[i] = (0, 0, 0, 0)
209#
210# rgb_alpha_image.putdata(pixel_data)
211# rgb_alpha_image.save(file_name, extension, compress_level=compress_level)
212
213# def create_raster(file_name: str, extension: str):
214# # Generating final raster image using concatenation.
215# image_data_dict = combine_image_data_by_x()
216# horizontally_concatenated_array = []
217#
218# # Concatenating horizontally each row.
219# for im_data_row in image_data_dict.values():
220# temp_array = im_data_row[0]
221# for im_index in range(1, len(im_data_row)):
222# temp_array = np.concatenate((temp_array, im_data_row[im_index]))
223# horizontally_concatenated_array.append(temp_array)
224#
225# # Concatenating vertically
226# final_resulting_array = horizontally_concatenated_array[0]
227# for row in range(1, len(horizontally_concatenated_array)):
228# final_resulting_array = np.concatenate((final_resulting_array,
229# horizontally_concatenated_array[row]), axis=1)
230# image_from_array(final_resulting_array, file_name, extension)