While using the new calculate_default_transform from Terracotta seemed to help with the performance, I encountered a but when dealing with scene that crosses the dateline separation

The problem happens here: https://github.com/cogeotiff/rio-tiler/blob/06504ec8321d185602a5a91bc6b5415bdb94915c/rio_tiler/utils.py#L232-L241

which will return negative width/height

Tile
x=17
y=21
z=6

vrt_transfrom
|-386.76, 0.00,-9392582.04|
| 0.00, 6.31, 6887893.49|
| 0.00, 0.00, 1.00|

vrt_width=-1619
vrt_height=-99186

dataset link: https://www.dropbox.com/s/jyskcxqnbq3240c/ask_cog.tif?dl=0

Not sure if it means something but here is the result from both rasterio and terracotta's calculate_default_transform for a random tile.

# rasterio
| 238.34, 0.00,-20037508.08|
| 0.00,-238.34, 13596067.67|
| 0.00, 0.00, 1.00|

# terracotta
|-386.54, 0.00,-20028306.76|
| 0.00, 6.31, 13596067.67|
| 0.00, 0.00, 1.00|

@dionhaefner, could you test terracotta with the dataset I shared, so I can confirm this is only happening in rio-tiler ?

https://github.com/DHI-GRAS/terracotta is a new lib inspired by rio-tiler made by @dionhaefner et all. This lib has plenty of nice tricks which could be useful in rio-tiler eg:

• different resampling method for up and down sampling: https://github.com/DHI-GRAS/terracotta/blob/7aa0751a105d1a804f0c5389e0abaab4399ffdee/terracotta/drivers/raster_base.py#L431-L436
• check tile coverage https://github.com/DHI-GRAS/terracotta/blob/7aa0751a105d1a804f0c5389e0abaab4399ffdee/terracotta/drivers/raster_base.py#L428-L429
• prevent interpolation artefact ? https://github.com/DHI-GRAS/terracotta/blob/7aa0751a105d1a804f0c5389e0abaab4399ffdee/terracotta/drivers/raster_base.py#L396-L411

I got this error: https://github.com/mapbox/rasterio/blob/master/rasterio/_warp.pyx#L907

and this code explains that it checks alpha band by GDAL color interpretation. https://github.com/mapbox/rasterio/blob/master/rasterio/_warp.pyx#L898

p.s: sorry for my bad english

Hello, i have a tiff file that i can open with rasterio.

But i get the following error with COGReader:

File "/usr/local/lib/python3.9/site-packages/rio_tiler/io/cogeo.py", line 101, in attrs_post_init self.nodata = self.nodata if self.nodata is not None else self.dataset.nodata AttributeError: '_GeneratorContextManager' object has no attribute 'nodata'

first mentioned in https://github.com/cogeotiff/rio-cogeo/issues/164

I'm seeing some weird Internal Tile fetching behavior. When working with a Web Mercator aligned file I would expect GDAL to only do 2 GET requests (header + tile) but It seems that we are fetching more than 1 tile.

$ tilebench profile https://rio-tiler-dev.s3.amazonaws.com/data/eu_webAligned_256pxWEBP.tif --tile 5-10-9 --add-kernels | jq
{
  "LIST": {
    "count": 0
  },
  "HEAD": {
    "count": 1
  },
  "GET": {
    "count": 3,
    "bytes": 17144,
    "ranges": [
      "0-16383",
      "205900-206279",
      "242386-242765"
    ]
  },
  "WarpKernels": [
    [
      "Src=512,2303,257x129",
      "Dst=0,0,256x128"
    ],
    [
      "Src=512,2431,257x129",
      "Dst=0,128,256x128"
    ]
  ],
  "Timing": 0.4853549003601074
}

☝️ we do 3 GET,

[
      "0-16383",  <-- GeoTIFF Header
      "205900-206279",  <--- Tile Data
      "242386-242765".  <--- Tile Data
]

The weird part is when we look at the WarpKernels:

[
    [
      "Src=512,2303,257x129",
      "Dst=0,0,256x128"
    ],
    [
      "Src=512,2431,257x129",
      "Dst=0,128,256x128"
    ]
]

I'm not sure how to read the Kernels but the 257x129 seems weird.

For the particular tile 5-10-9, here are the WarpedVRT info

{
    'add_alpha': True, 
    'resampling': <Resampling.nearest: 0>, 
    'crs': CRS.from_epsg(3857), 
    'transform': Affine(2445.984905125635, 0.0, -7514065.6285459455, 0.0, -2445.984905125635, 8766409.89997027), 
    'width': 512, 
    'height': 512
}
Window(col_off=0, row_off=0, width=512, height=512)
(3, 256, 256)

The VRT is 512x512 for a 256x256 tile, which is what we expect because the zoom 5 correspond to the first overview.

The input and output size are ok, so maybe we need to check if the transform is aligned with the internal tiles. https://github.com/cogeotiff/rio-tiler/blob/4bd9be5db98827e40d4222627c3af59e9063bd95/rio_tiler/utils.py#L166-L182

This is a major change

While working on #95 and https://github.com/RemotePixel/amazonlinux/issues/10 I realized that maybe we were using WarpedVRT in an unoptimized way, when (comparing warped kernel + GET requests)

This goes back to https://github.com/mapbox/rasterio/issues/1373 and precisely https://github.com/mapbox/rasterio/issues/1373#issuecomment-398487266

When looking at the gdalwarp command I realized the gdal use specific overview level to create the output dataset. This PR tries to implement such logic.

The first tests show that GET requests and output results are similar to gdalwarp outputs.

~@sgillies I'll need your help on this one (and I believe this will also fix the performance downgrade I saw in https://github.com/RemotePixel/amazonlinux/issues/10).~ (Edit the performance downgrade was due to SQLITE3 and centos) I'm not sure if this could also go directly to rasterio code base. (I know you have ton of stuff on your side so no stress, I'll continue to run some tests on my side).

tagging @dionhaefner @mojodna because it might be of your interest.

in https://github.com/cogeotiff/rio-tiler/blob/2bf9f144958ae226af4db338487e8a34159550be/rio_tiler/utils.py#L345-L353 we could end up (actual behavior) with using different resampling for mask and data which will result in bad masking.

Reopening issue #53 because this has not been fixed and is an issue with rio-tiler.

Running the following:

from rio_tiler import sentinel2
import os
os.environ["CURL_CA_BUNDLE"] = "/etc/ssl/certs/ca-certificates.crt"
os.environ["AWS_REQUEST_PAYER"] = "requester"
tile, mask = sentinel2.tile("S2B_tile_20180813_23MQQ_0", 781, 1052, 11, (4, 3, 2), tilesize=256)

gives the error:

RasterioIOError: '/vsis3/sentinel-s2-l1c/tiles/23/M/QQ/2018/8/13/0/B4.jp2' does not exist in the file system, and is not recognized as a supported dataset name.

running gdalinfo on this path gives the error:

ERROR 4: /vsis3/sentinel-s2-l1c/tiles/23/M/QQ/2018/8/13/0/B4.jp2: No such file or directory
gdalinfo failed - unable to open '/vsis3/sentinel-s2-l1c/tiles/23/M/QQ/2018/8/13/0/B4.jp2'.

This path is being generated by rio-tiler. This is only happening on Sentinel; Landsat works fine.

rio-tiler version: rio-tiler==2.0.0rc4

Recreate:

from rio_tiler.io import COGReader
import matplotlib.pyplot as plt


feature = ... # huge multipolygon I'll attach separately

with COGReader('https://rio-tiler-dev.s3.amazonaws.com/data/eu_webAligned_256pxWEBP.tif') as cog:
    data, mask = cog.feature(feature['geometry'])

plt.imshow(mask)
plt.show()

Tried the .buffer(0) trick as well with no change to the result. Geometry returns True for is_valid

example (EPSG:32601)

Upper Left  (365474, 7696071) 
Lower Left  (365474, 7591951) 
Upper Right (477094, 7696071)
Lower Right (477094, 7591951)

rasterio. warp.transform_bounds return flattened bounds.

transform_bounds('EPSG:32601', 'EPSG:4326', 365474, 7591951, 477094, 7696071)
-> (-179.90689695793938, 68.40816857051695, 179.96969316611757, 69.37306968939305)

correct bounds (179.5819541972332, 68.40816857051695, -177.55854062377162, 69.37306968939305)

It seems to need checking transformed points is crossing the antimeridian. For example, like this.

def is_clockwise(xs, ys):
    if len(xs) != len(ys):
        raise TransformError("xs and ys arrays must be the same length") 
    
    area = 0
    if len(xs) > 2:
        if xs[0] == xs[-1] and ys[0] == ys[-1]:
            length = len(xs) - 1
        else:
            length = len(xs)
        
        for i in range(length - 1):
            area = area + xs[i] * ys[i+1] - ys[i] * xs[i+1]
        area = area + xs[length - 1] * ys[0] - ys[length - 1] * xs[0]
            
    return area >= 0

def transform_bounds(
        src_crs,
        dst_crs,
        left,
        bottom,
        right,
        top,
        densify_pts=21):
    if densify_pts < 0:
        raise ValueError('densify parameter must be >= 0')

    if src_crs == dst_crs:
        return [left, bottom, right, top]
    
    in_xs = []
    in_ys = []

    if densify_pts > 0:
        densify_factor = 1.0 / float(densify_pts + 1)
        
        # left_bottom to right_bottom
        in_xs.extend(
            left + np.arange(1, densify_pts + 1, dtype=np.float64) *
            ((right - left) * densify_factor)
        )
        in_ys.extend([bottom] * densify_pts)
        # right_bottom to right_top
        in_xs.extend([right] * (densify_pts + 2))
        in_ys.extend(
            bottom + np.arange(0, densify_pts + 2, dtype=np.float64) *
            ((top - bottom) * densify_factor)
        )
        # right_top to left_top
        in_xs.extend(
            right + np.arange(1, densify_pts + 1, dtype=np.float64) *
            ((left - right) * densify_factor)
        )
        in_ys.extend([top] * densify_pts)
        # left_top to left_bottom
        in_xs.extend([left] * (densify_pts + 2))
        in_ys.extend(
            top + np.arange(0, densify_pts + 2, dtype=np.float64) *
            ((bottom - top) * densify_factor)
        )        

    else:
        in_xs = np.array(left, right, right, left,left)
        in_ys = np.array(bottom, bottom, top, top, bottom)

    xs, ys = transform(src_crs, dst_crs, in_xs, in_ys)
    if dst_crs == 'EPSG:4326' and not is_clockwise(xs, ys):
        xs = [x + 360 if x < 0 else x for x in xs]
        return (min(xs), min(ys), max(xs) - 360, max(ys))
    
    return (min(xs), min(ys), max(xs), max(ys))

cogeo.bounds, cogeo.get_zooms, reader.part , etc., a lot of methods to create tile are affected. Couldn't we improve?

As proposed in #56 this PR aims to reduce the edge effects associated with tiles by introducing a tile_edges_padding parameter. Additional data is retrieved around the tiles edges during the tile generation process which helps to reduce sharp edges that occur due to resampling.

@vincentsarago did suggest

and then overwrite it to 0 if we don't need resampling (if the internal tiles match the requested bounds)

I'm not really sure how we'd know about the internal tiles so happy to take any advice on this one :)

Let me know if you need anything else on this.

Let's say you have an Non-geo image (as COG already) and some external GCP (e.g from a Georeferencing UI), it would be nice to pass those GCPS to rio_tiler.io.Reader directly to allow virtual georeferencing

https://github.com/cogeotiff/rio-tiler/blob/93cf68ed9d4ba62c9fb51a11f601c1a141324ebb/rio_tiler/io/rasterio.py#L98-L110

Let's make it clear what user should write to avoid issues like #546 #545

ref: https://docs.github.com/en/communities/using-templates-to-encourage-useful-issues-and-pull-requests

Here is a simple template we have at devseed

<!-- Please search existing issues to avoid creating duplicates. -->


#### Problem description

[this should explain **why** the current behavior is a problem and why the expected output is a better solution.]

#### Expected Output


#### Environment Information

Os, Python version, ...

There is following set of files:

• masked.tif (with PER_DATASET internal mask)
• data.vrt built based on masked.tif
• data.vrt.ovr (external overviews)

Let's try to render a tile out of it:

In [1]: from rio_tiler.io import COGReader

In [2]: with COGReader("/tmp/test/data.vrt") as cog:
   ...:     tile = cog.tile(2323, 3503, 13, tilesize=256)
   ...: 
GDAL: GDALClose(/tmp/test/masked.tif, this=0x3699940)
GDAL: GDALClose(/tmp/test/data.vrt, this=0x3640b10)

In [3]: with open("/tmp/test/tile-ok.png", "wb") as dst:
   ...:     dst.write(tile.render())

Expected result (screenshot from QGIS):

Actual result:

Everything is fine.

Now let's try to render tile on a lower zoom level, in this case GDAL will use existing overviews.

In [1]: from rio_tiler.io import COGReader

In [2]: with COGReader("/tmp/test/data.vrt") as cog:
   ...:     tile = cog.tile(72, 109, 8, tilesize=256)
   ...: 
GDAL: GDALClose(/tmp/test/masked.tif, this=0x2f1f7e0)
GDAL: GDALClose(/tmp/test/data.vrt, this=0x20b6740)
GDAL: GDALClose(/tmp/test/data.vrt.ovr, this=0x2e51cc0)

In [3]: with open("/tmp/test/tile-nok.png", "wb") as dst:
   ...:     dst.write(tile.render())
   ...: 

Expected result:

Actual result:

As you can see the mask has not been taken into account.

Files that are being used: example.zip

By design rio-tiler assume you are working within the native zoom range of your data (raw data = Max res, max res/number of overviews ~= minzoom).

When reading a file, we always use a VRT to align the data to the bounds/size we want to read. This is all defined in https://github.com/cogeotiff/rio-tiler/blob/ccf6aa002fc0c3525c9889cffbb9ca26689ba0b9/rio_tiler/utils.py#L172-L188

(this part was copie from https://github.com/mapbox/rasterio/issues/1373#issuecomment-398261502)

if you try to read a really low version of your data, we will first create a huge VRT and then ask to read a decimated part of it (this is to make sure we read the overview).

To be honest, for now I'm not sure there is a solution for this.

Using img.render() with an ImageData.data array with too many bands fails with CPLE_NotSupportedError: PNG driver doesn't support 5 bands. Must be 1 (grey), 2 (grey+alpha), 3 (rgb) or 4 (rgba) bands.

Example:

src_path = "https://naipblobs.blob.core.windows.net/naip/v002/al/2019/al_60cm_2019/30087/m_3008701_ne_16_060_20191115.tif"
with COGReader(src_path) as cog:
    img = cog.preview()    
    print(img.data.shape)  # >>> (4, 1024, 892)

img.render()

Reading all bands and later visualizing a band subset (without mutating the array) seems like a common usecase. Would propose to add an indexes parameter to .render(), so that specific bands can be selected for the vizualization.

Happy to make a PR, but wanted to check first if there is a better solution or this is not intended etc.