force overview dataset for WarpedVRT #132
Closed
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
vincentsarago
commented
Nov 22, 2019
| options = {"overview_level": overview_level} if overview_level is not None else {} | ||
| with rasterio.open(src_dst.name, **options) as ovr_dst: | ||
| w, s, e, n = bounds | ||
| vrt_transform = transform.from_bounds(w, s, e, n, tilesize, tilesize) |
There was a problem hiding this comment.
This is the main change.
Instead of relying on WarpedVRT, we force rasterio to use a specific overview level (or the highest resolution)
Note, we cannot pass overview_level=None to rasterio because it somehow get converted to overview_level=0
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
This comment was marked as resolved.
This comment was marked as resolved.
Sorry, something went wrong.
|
Notes: Old Rio-TilerNew Rio-TilerGDAL_DISABLE_READDIR_ON_OPENNote:
import sys
import logging
import mercantile
import rasterio
from rasterio import transform
from rio_tiler.utils import array_to_image, _tile_read
logging.basicConfig(stream=sys.stderr, level=10)
def _rio_tiler_vrt(resampling, bounds, tilesize=256):
with rasterio.open(input) as src_dst:
tile, _ = _tile_read(
src_dst,
bounds,
tilesize=tilesize,
resampling_method=resampling,
tile_edge_padding=0,
)
dst_crs = {'init': 'EPSG:3857'}
w, s, e, n = bounds
dst_transform = transform.from_bounds(w, s, e, n, tilesize, tilesize)
with open(f"{tile_z}-{tile_x}-{tile_y}_{resampling}.tif", "wb") as f:
f.write(array_to_image(
tile,
dtype=tile.dtype,
img_format="GTiff",
crs=dst_crs,
transform=dst_transform
))
if __name__ == '__main__':
input = "https://s3-us-west-2.amazonaws.com/remotepixel-pub/data/Int16_nodata9999.tif"
tile_z = 12
tile_x = 676
tile_y = 1619
tile = mercantile.Tile(x=tile_x, y=tile_y, z=tile_z)
tile_bounds = mercantile.xy_bounds(tile)
_rio_tiler_vrt("bilinear", tile_bounds, tilesize=256) |
|
Note2 (fetching overview) Old rio-tilerNew Rio-tilerGDALAgain, both new and old fetch the same amount of data but the new rio-tiler gives the exact same result as the gdalwarp command! |
|
Note: started a quick benchmarking suite over https://github.com/vincentsarago/rio-tiler-bench For now it shows real good performances with this new update |
|
Note: The tests are failing because we are trying to re-open a WarpedVRT dataset (sentinel-1). Not sure what's the best path but maybe I could spend some time to see if rasterio/rasterio#1504 can be a solution. with rasterio.open(src) as src_dst:
ovr = get_overview_level(...)
with rasterio.open(src_dst.name, overview_level=ovr) as ovr_dst:
.... we could do with rasterio.open(src) as src_dst:
ovr = get_overview_level(...)
ovr_dst = src_dst.overview(ovr) if ovr else src_dst
.... |
|
Is this only an issue for rasters with an internal mask? |
@dionhaefner no it's not, this is more a general fix because I think the way we were building the WarpedVRT before is not optimal! The main motivation was to fix performances saw with GDAL3 but also to try to get the same results than with gdalwarp commands |
This comment has been minimized.
This comment has been minimized.
|
Benchmark results on Terracotta: (GDAL 2) Current masterPassing
|
dionhaefner
reviewed
Nov 26, 2019
|
Well I'm going to close this. Rasterio WarpedVRT is doing the right thing and fetch the overview, the problem I was seeing (difference in kernel) due to the |
|
🤔 this is not over in fact, I was able to confirm that forcing the process on the overview was giving the same result as with non-force overview for -> https://gist.github.com/vincentsarago/9951d1aab15359ee826b3b4c2cebeeab @sgillies This might be a pure resampling problem and not really a WarpedVRT thing. For now I'm going to remove the padding option (which was causing most of the problem) and document the use of SOURCE_EXTRA |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This is a major change
While working on #95 and RemotePixel/amazonlinux#10 I realized that maybe we were using WarpedVRT in an
unoptimizedway, when (comparing warped kernel + GET requests)This goes back to rasterio/rasterio#1373 and precisely rasterio/rasterio#1373 (comment)
When looking at the
gdalwarpcommand 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 RemotePixel/amazonlinux#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.