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