Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
arr = numpy.transpose(arr, [2, 0, 1])
tbounds = mercantile.xy_bounds(leftTile)
data, mask = tile_read(
"cogeo.tif", tbounds, 256, resampling_method="nearest"
)
numpy.testing.assert_array_equal(data, arr)
# Bottom right tile
mime_type, tile = cog.get_tile(4, 3, 0)
tile_length = 256 * 256 * 3
t = struct.unpack_from("{}b".format(tile_length), tile)
arr = numpy.array(t).reshape(256, 256, 3).astype(numpy.uint8)
arr = numpy.transpose(arr, [2, 0, 1])
tbounds = mercantile.xy_bounds(rightTile)
data, mask = tile_read(
"cogeo.tif", tbounds, 256, resampling_method="nearest"
)
numpy.testing.assert_array_equal(data, arr)
# Low resolution (overview 1)
# Top Left tile
# NOTE: overview internal tile size is 128px
# We need to stack two internal tiles to compare with
# the 256px mercator tile fetched by rio-tiler
# ref: https://github.com/cogeotiff/rio-cogeo/issues/60
mime_type, tile = cog.get_tile(1, 0, 1)
tile_length = 128 * 128 * 3
t = struct.unpack_from("{}b".format(tile_length), tile)
arr1 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
arr1 = numpy.transpose(arr1, [2, 0, 1])
def read_tile(self, z, x, y):
"""Read raster tile data and mask."""
mercator_tile = mercantile.Tile(x=x, y=y, z=z)
tile_bounds = mercantile.xy_bounds(mercator_tile)
data, mask = tile_read(
self.path,
tile_bounds,
self.tiles_size,
indexes=self.indexes,
nodata=self.nodata,
)
data = (data[0] + data[1]) / 2
return data.astype(numpy.uint8), mask
im = Image.fromarray(img_arr)
x = 0
y = 0
im_path = join(self.tmp_dir, '{}/{}/{}.png'.format(zoom, x, y))
make_dir(im_path, use_dirname=True)
im.save(im_path)
tile_schema = join(self.tmp_dir, '{z}/{x}/{y}.png')
# Get NW and SE corner of central half of tile.
nw_bounds = mercantile.xy_bounds(0, 0, zoom)
nw_merc_y = nw_bounds.top - (nw_bounds.top - nw_bounds.bottom) / 4
nw_merc_x = nw_bounds.left + (nw_bounds.right - nw_bounds.left) / 4
nw_lng, nw_lat = merc2lnglat(nw_merc_x, nw_merc_y)
se_bounds = mercantile.xy_bounds(0, 0, zoom)
se_merc_y = se_bounds.bottom + (se_bounds.top - se_bounds.bottom) / 4
se_merc_x = se_bounds.right - (se_bounds.right - se_bounds.left) / 4
se_lng, se_lat = merc2lnglat(se_merc_x, se_merc_y)
# min_lat, min_lng, max_lat, max_lng = bounds
bounds = [se_lat, nw_lng, nw_lat, se_lng]
output_uri = join(self.tmp_dir, 'output.tif')
_zxy2geotiff(tile_schema, zoom, bounds, output_uri)
with rasterio.open(output_uri) as dataset:
tiff_arr = dataset.read()
self.assertEqual(tiff_arr.shape, (3, 128, 128))
exp_arr = np.transpose(img_arr, (2, 0, 1))[:, 64:-64, 64:-64]
np.testing.assert_array_equal(tiff_arr, exp_arr)
# ref: https://github.com/cogeotiff/rio-cogeo/issues/60
mime_type, tile = cog.get_tile(1, 0, 1)
tile_length = 128 * 128 * 3
t = struct.unpack_from("{}b".format(tile_length), tile)
arr1 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
arr1 = numpy.transpose(arr1, [2, 0, 1])
mime_type, tile = cog.get_tile(2, 0, 1)
tile_length = 128 * 128 * 3
t = struct.unpack_from("{}b".format(tile_length), tile)
arr2 = numpy.array(t).reshape(128, 128, 3).astype(numpy.uint8)
arr2 = numpy.transpose(arr2, [2, 0, 1])
arr = numpy.dstack((arr1, arr2))
lowTile = mercantile.Tile(118594, 60034, 17)
tbounds = mercantile.xy_bounds(lowTile)
data, mask = tile_read(
"cogeo.tif", tbounds, 256, resampling_method="nearest"
)
data = data[:, 128:, :]
numpy.testing.assert_array_equal(data, arr)
blocks = list(set(out_dst.block_shapes))
assert len(blocks) == 1
ts = blocks[0][0]
assert not out_dst.width % ts
assert not out_dst.height % ts
max_zoom = get_max_zoom(out_dst)
bounds = list(
transform_bounds(
*[src_dst.crs, "epsg:4326"] + list(src_dst.bounds),
densify_pts=21
)
)
leftTile = mercantile.tile(bounds[0], bounds[3], max_zoom)
tbounds = mercantile.xy_bounds(leftTile)
west, north = tbounds.left, tbounds.top
assert out_dst.transform.xoff == west
assert out_dst.transform.yoff == north
rightTile = mercantile.tile(bounds[2], bounds[1], max_zoom)
tbounds = mercantile.xy_bounds(rightTile)
east, south = tbounds.right, tbounds.bottom
lrx = round(
out_dst.transform.xoff + out_dst.transform.a * out_dst.width, 6
)
lry = round(
out_dst.transform.yoff + out_dst.transform.e * out_dst.height, 6
)
assert lrx == round(east, 6)
assert lry == round(south, 6)
def _expand_bounds(self, bounds):
if bounds is None:
return bounds
min_tile_x, min_tile_y, max_tile_x, max_tile_y = self._tile_coords(bounds)
ul = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=min_tile_x, y=max_tile_y)))
lr = box(*mercantile.xy_bounds(mercantile.Tile(z=self.zoom_level, x=max_tile_x, y=min_tile_y)))
return ul.union(lr).bounds
def merc2pixel(tile_x, tile_y, zoom, merc_x, merc_y, tile_sz=256):
"""Convert Web Mercator point to pixel coordinates.
This is within the coordinate frame of a single ZXY tile.
Args:
tile_x: (int) x coordinate of ZXY tile
tile_y: (int) y coordinate of ZXY tile
zoom: (int) zoom level of ZXY tile
merc_x: (float) Web Mercator x axis of point
merc_y: (float) Web Mercator y axis of point
tile_sz: (int) size of ZXY tile
"""
tile_merc_bounds = mercantile.xy_bounds(tile_x, tile_y, zoom)
pix_y = int(
round(tile_sz * ((tile_merc_bounds.top - merc_y) /
(tile_merc_bounds.top - tile_merc_bounds.bottom))))
pix_x = int(
round(tile_sz * ((merc_x - tile_merc_bounds.left) /
(tile_merc_bounds.right - tile_merc_bounds.left))))
return (pix_x, pix_y)
"""Render a tile into Web Mercator.
Arguments:
tile {mercantile.Tile} -- Tile to render.
sources {list} -- Sources to render from.
Keyword Arguments:
transformation {Transformation} -- Transformation to apply. (default: {None})
format {function} -- Output format. (default: {None})
scale {int} -- Output scale factor. (default: {1})
expand {bool} -- Whether to expand single-band, paletted sources to RGBA. (default: {True})
Returns:
(dict, bytes) -- Tuple of HTTP headers (dict) and bytes.
"""
bounds = Bounds(mercantile.xy_bounds(tile), WEB_MERCATOR_CRS)
shape = tuple(map(int, Affine.scale(scale) * TILE_SHAPE))
return render(
bounds,
shape,
WEB_MERCATOR_CRS,
sources=sources,
format=format,
transformation=transformation,
expand=expand,
)
if band not in LANDSAT_BANDS:
raise InvalidBandName("{} is not a valid Landsat band name".format(band))
scene_params = _landsat_parse_scene_id(sceneid)
meta_data = _landsat_get_mtl(sceneid).get("L1_METADATA_FILE")
landsat_address = "{}/{}".format(LANDSAT_BUCKET, scene_params["key"])
wgs_bounds = toa_utils._get_bounds_from_metadata(meta_data["PRODUCT_METADATA"])
if not utils.tile_exists(wgs_bounds, tile_z, tile_x, tile_y):
raise TileOutsideBounds(
"Tile {}/{}/{} is outside image bounds".format(tile_z, tile_x, tile_y)
)
mercator_tile = mercantile.Tile(x=tile_x, y=tile_y, z=tile_z)
tile_bounds = mercantile.xy_bounds(mercator_tile)
def _tiler(band):
address = "{}_B{}.TIF".format(landsat_address, band)
if band == "QA":
nodata = 1
else:
nodata = 0
return utils.tile_read(
address, bounds=tile_bounds, tilesize=tilesize, nodata=nodata, **kwargs
)
with futures.ThreadPoolExecutor(max_workers=MAX_THREADS) as executor:
data, masks = zip(*list(executor.map(_tiler, bands)))
data = np.concatenate(data)
mask = np.all(masks, axis=0).astype(np.uint8) * 255
"""Burn tile with features.
Args:
tile: the mercantile tile to burn.
features: the geojson features to burn.
size: the size of burned image.
Returns:
image: rasterized file of size with features burned.
"""
# the value you want in the output raster where a shape exists
burnval = 1
shapes = ((geometry, burnval) for feature in features for geometry in feature_to_mercator(feature))
bounds = mercantile.xy_bounds(tile)
transform = from_bounds(*bounds, size, size)
return rasterize(shapes, out_shape=(size, size), transform=transform)