1

I am trying to use rasterio to read cloud-optimized geotiff stored in a gc bucket.

With method 1, GDAL makes use of the overviews and read is quick. Unfortunately, if the extent is not fully included in the original raster, data is "strechted", like if extent was cropped to the projected raster extent, then stretch to fill the out_shape array, forgetting the no_data value.

image

With method 2, read is extremely long as it looks like overview reads are no triggered (several hundreds of requests).

What is the right way to read cogeo in the cloud ?

# tile
z = 8
x = 126
y = 202
b = 15

# tile properties
wm = web_mercator(z)
extent = wm.tile_extent_xy(a=x, b=y)
transform = wm.tile_transform(a=x, b=y)

# method 1 - quick
# - triggers the use of overviews
# - incorrect data at borders - data is strechted if extent not fully included in raster
profile = {
    'crs': 'epsg:3857',
    'resampling': rasterio.enums.Resampling.nearest
}

with rasterio.Env(**gdal_env):
    with rasterio.open(raster) as src:
        with rasterio.vrt.WarpedVRT(src, **profile) as vrt:
            copy(vrt, 'test1.vrt', driver='vrt')
            window = vrt.window(*extent)
            data = vrt.read(b,
                            window=window,
                            out_shape=(256, 256))

# method 2 - slow
# - get all data at base resolution
profile = {
    'crs': 'epsg:3857',
    'transform': transform,
    'width': 256,
    'height': 256,
    'resampling': rasterio.enums.Resampling.nearest,
    'src_nodata': None
}

with rasterio.Env(**gdal_env):
    with rasterio.open(raster) as src:
        with rasterio.vrt.WarpedVRT(src, **profile) as vrt:
            copy(vrt, 'test2.vrt', driver='vrt')
            data = vrt.read(b)

Here's test1.xml :

<VRTDataset rasterXSize="11251" rasterYSize="10421" subClass="VRTWarpedDataset">
  <SRS>PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</SRS>
  <GeoTransform>  1.1478038722865377e+07,  2.0073715778544621e+01,  0.0000000000000000e+00,  2.6634411708228581e+05,  0.0000000000000000e+00, -2.0073715778544621e+01</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1" subClass="VRTWarpedRasterBand">
    <ColorInterp>Gray</ColorInterp>
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="2" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="3" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="4" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="5" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="6" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="7" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="8" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
    <ResampleAlg>NearestNeighbour</ResampleAlg>
    <WorkingDataType>Int16</WorkingDataType>
    <Option name="UNIFIED_SRC_NODATA">YES</Option>
    <Option name="INIT_DEST">NO_DATA</Option>
    <SourceDataset relativeToVRT="0">/vsigs/project/P125R59_L5_1993.tif</SourceDataset>
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGeoTransform>289720,20,0,264380,0,-20</SrcGeoTransform>
            <SrcInvGeoTransform>-14486,0.0500000000000000028,0,13219,0,-0.0500000000000000028</SrcInvGeoTransform>
            <DstGeoTransform>11478038.7228653766,20.0737157785446207,0,266344.117082285811,0,-20.0737157785446207</DstGeoTransform>
            <DstInvGeoTransform>-571794.422591827461,0.0498163873112535266,0,13268.301694645008,0,-0.0498163873112535266</DstInvGeoTransform>
            <ReprojectTransformer>
              <ReprojectionTransformer>
                <SourceSRS>PROJCS["WGS 84 / UTM zone 48N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32648"]]</SourceSRS>
                <TargetSRS>PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</TargetSRS>
              </ReprojectionTransformer>
            </ReprojectTransformer>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1" />
      <BandMapping src="2" dst="2" />
      <BandMapping src="3" dst="3" />
      <BandMapping src="4" dst="4" />
      <BandMapping src="5" dst="5" />
      <BandMapping src="6" dst="6" />
      <BandMapping src="7" dst="7" />
      <BandMapping src="8" dst="8" />
    </BandList>
  </GDALWarpOptions>
</VRTDataset>

And test2.vrt :

<VRTDataset rasterXSize="256" rasterYSize="256" subClass="VRTWarpedDataset">
  <SRS>PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</SRS>
  <GeoTransform>  1.1584184510675032e+07,  6.1149622628141003e+02,  0.0000000000000000e+00,  3.1308606785608083e+05,  0.0000000000000000e+00, -6.1149622628141003e+02</GeoTransform>
  <VRTRasterBand dataType="Int16" band="1" subClass="VRTWarpedRasterBand">
    <ColorInterp>Gray</ColorInterp>
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="2" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="3" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="4" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="5" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="6" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="7" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <VRTRasterBand dataType="Int16" band="8" subClass="VRTWarpedRasterBand">
  </VRTRasterBand>
  <BlockXSize>256</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <WarpMemoryLimit>6.71089e+07</WarpMemoryLimit>
    <ResampleAlg>NearestNeighbour</ResampleAlg>
    <WorkingDataType>Int16</WorkingDataType>
    <Option name="UNIFIED_SRC_NODATA">YES</Option>
    <Option name="INIT_DEST">NO_DATA</Option>
    <SourceDataset relativeToVRT="0">/vsigs/project/P125R59_L5_1993.tif</SourceDataset>
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGeoTransform>289720,20,0,264380,0,-20</SrcGeoTransform>
            <SrcInvGeoTransform>-14486,0.0500000000000000028,0,13219,0,-0.0500000000000000028</SrcInvGeoTransform>
            <DstGeoTransform>11584184.5106750317,611.496226281410031,0,313086.06785608083,0,-611.496226281410031</DstGeoTransform>
            <DstInvGeoTransform>-18944,0.0016353330683349153,0,511.999999999998181,0,-0.0016353330683349153</DstInvGeoTransform>
            <ReprojectTransformer>
              <ReprojectionTransformer>
                <SourceSRS>PROJCS["WGS 84 / UTM zone 48N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32648"]]</SourceSRS>
                <TargetSRS>PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]</TargetSRS>
              </ReprojectionTransformer>
            </ReprojectTransformer>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1" />
      <BandMapping src="2" dst="2" />
      <BandMapping src="3" dst="3" />
      <BandMapping src="4" dst="4" />
      <BandMapping src="5" dst="5" />
      <BandMapping src="6" dst="6" />
      <BandMapping src="7" dst="7" />
      <BandMapping src="8" dst="8" />
    </BandList>
  </GDALWarpOptions>
</VRTDataset>
Improve this question

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.