Use QuadBin - Githubissues

francois-baptiste commented 1 year ago

Here is a fixture tiff https://github.com/CartoDB/raster-loader/blob/661adc9a562792a60bf59edf3a542d2af7aff242/raster_loader/tests/fixtures/world_3857.tif

DEPRECATED: see https://github.com/CartoDB/raster-loader/issues/45#issuecomment-1353288934

brendancol commented 1 year ago

@francois-baptiste Nice...XYZ tiles would also involve a resampling method. Would you imagine a table like: | X | Y | Z | landcat_uint8 |

|0 | 0 | 0 | Bytes()|

francois-baptiste commented 1 year ago

Yes exactly.

In a first time, if we support only EPSG Code: 3857 (webmercator) we do not need reprojection. There is a direct match between:

each pixel and a (X,Y,Z) tuple
each pixel blocs and a (X,Y,Z) tuple

brendancol commented 1 year ago

We will assume WMTS tile structure

brendancol commented 1 year ago

We just need to store X,Y,Z columns that can later be converted into Carto Quadbin

brendancol commented 1 year ago

we will assume the data EXACTLY matches the resolution of a standard web mercator level

francois-baptiste commented 1 year ago

Thanks to @jatorre I have generated a tiff feature that should match exactly the resolution of a standard web mercator level

https://github.com/CartoDB/raster-loader/blob/61ccde5dee6d2b01b102960e1a2b3753ab415e9e/raster_loader/tests/fixtures/world_zoom_0-2.cog.tif

Command I used to generate it (be sure to be on the latest version of gdal):

gdalwarp -of COG -co TILING_SCHEME=GoogleMapsCompatible -co COMPRESS=DEFLATE world_zoom_0-2.mbtiles world_zoom_0-2.cog.tif`

gdalinfo /home/user/Downloads/world_zoom_0-2.cog.tif
Driver: GTiff/GeoTIFF
Files: /home/user/Downloads/world_zoom_0-2.cog.tif
Size is 1024, 1024
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-20037508.342789243906736,20037508.342789243906736)
Pixel Size = (39135.758482010242005,-39135.758482010242005)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  LAYOUT=COG
Tiling Scheme:
  NAME=GoogleMapsCompatible
  ZOOM_LEVEL=2
Corner Coordinates:
Upper Left  (-20037508.343,20037508.343) (180d 0' 0.00"W, 85d 3' 4.06"N)
Lower Left  (-20037508.343,-20037508.343) (180d 0' 0.00"W, 85d 3' 4.06"S)
Upper Right (20037508.343,20037508.343) (180d 0' 0.00"E, 85d 3' 4.06"N)
Lower Right (20037508.343,-20037508.343) (180d 0' 0.00"E, 85d 3' 4.06"S)
Center      (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 512x512, 256x256
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 512x512, 256x256
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 512x512, 256x256
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 512x512, 256x256
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 512x512, 256x256
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 512x512, 256x256
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 512x512, 256x256

francois-baptiste commented 1 year ago

Her is a POC for world_zoom_0-2.cog.tif Quadbin aggregation level 2 Quadbin pixel level 5

import pandas as pd
import rasterio
from google.cloud import bigquery

client = bigquery.Client()

with rasterio.open('/home/user/Downloads/world_zoom_0-2.cog.tif') as dataset:

    toto=dataset.read(1) #read raster band 1

    mydf = pd.DataFrame([(2, # TODO Rasterio should helps us to find a way to infer the Z level 
                          *ij, # X and Y level TOCHECK maybe I swap those or inverse Y
                          window.height,window.width, # maybe in a near future we can skip these columns and always assume 256x256
                          dataset.read(1, window=window).tobytes())
    for ij, window in dataset.block_windows()],columns=['Z','X','Y','block_height','block_width','value'])

print(mydf) #print dataframe

job=client.load_table_from_dataframe(mydf, 'yourtable.yourdataset.world_zoom_2')
job.result()

returns:

     Z  X  Y  block_height  block_width                                              value
0   2  0  0           256          256  b'\xd4\xd3\xd3\xd3\xd2\xd3\xd5\xd8\xd8\xd8\xd7...
1   2  0  1           256          256  b'\xd6\xd6\xd8\xd7\xd7\xd5\xd4\xd4\xd5\xd5\xd5...
2   2  0  2           256          256  b"\xcb\xcb\xcc\xcb\xc7\xcb\xd2\xd3\xd3\xc7\xbe...
3   2  0  3           256          256  b'\xcd\xce\xce\xce\xcc\xcb\xcd\xcc\xcc\xcd\xcd...
4   2  1  0           256          256  b"\xcf\xc4\xc4\xcc\xd5\xd9\xdb\xdc\xdb\xdb\xd9...
5   2  1  1           256          256  b'\xe1\xe1\xe1\xe4\xe1\xe1\xe1\xe4\xe0\xdf\xda...
6   2  1  2           256          256  b'\xce\xce\xce\xce\xd1\xd7\xd8\xd8\xd8\xd8\xda...
7   2  1  3           256          256  b'\xba\xb5\xba\xbe\xb8\xb7\xb3\xb1\xb3\xb8\xba...
8   2  2  0           256          256  b"\xbb\xb9\xbc\xbd\xbd\xb7\xbd\xbd\xbc\xbc\xb5...
9   2  2  1           256          256  b'\xc8\xca\xc9\xc8\xce\xd2\xcd\xd2\xd3\xd3\xd4...
10  2  2  2           256          256  b"\xbd\xbd\xbd\xbc\xbb\xbc\xbb\xb6\xbc\xbc\xbd...
11  2  2  3           256          256  b'\xc6\xbd\xbb\xbc\xbd\xbe\xbe\xbe\xbe\xbe\xbd...
12  2  3  0           256          256  b'\xcb\xc8\xc9\xcb\xc2\xc1\xc9\xc3\xce\xd1\xcf...
13  2  3  1           256          256  b'\xbd\xbd\xbd\xbd\xbd\xbd\xbe\xbe\xbd\xbd\xbd...
14  2  3  2           256          256  b'\xbb\xca\xc9\xbd\xc0\xbe\xc0\xca\xcb\xca\xcc...
15  2  3  3           256          256  b'\xd3\xd2\xe1\xfa\xf9\xf7\xf6\xee\xf8\xf8\xfa...

francois-baptiste commented 1 year ago

Right now the neatest approach I've found to infer quadbin cells is to use rio_cogeo to get max zoom level. Some error handling should be added though...

import pandas as pd
import quadbin
import pyproj
import rasterio
import rio_cogeo

raster_info = rio_cogeo.cog_info('world_zoom_0-2.cog.tif').dict()
max_zoom = raster_info["GEO"]["MaxZoom"]

dst_crs = pyproj.CRS.from_epsg(4326)

with rasterio.open('world_zoom_0-2.cog.tif') as dataset:

    src_crs = pyproj.CRS.from_wkt(dataset.crs.to_wkt())
    transformer = pyproj.Transformer.from_crs(src_crs, dst_crs) #compute lat and lon

    mydf = pd.DataFrame([(
        quadbin.point_to_cell(
            *reversed(
                transformer.transform(
                    *rasterio.transform.xy(
                        dataset.transform,
                        window.row_off+window.height*.5,
                        window.col_off+window.width*.5))),
            max_zoom),
        window.height,
        window.width,
        *idx,
        dataset.read(1, window=window).tobytes()
        )
        for idx, window
        in dataset.block_windows()],
        columns=['quadbin','block_height','block_width','block_height_idx','block_width_idx','value_uint8'])

francois-baptiste commented 1 year ago

Here is the definition of GoogleMapsCompatible Tiling scheme

https://github.com/OSGeo/gdal/blob/0a1eb9bd7aad7994370a89a58e3ade744475be12/doc/source/drivers/raster/gpkg.rst#tiling-schemes

francois-baptiste commented 1 year ago

@brendancol I have started a draft pull request for this issue. It would be great if you could find some time to review it once I've finalized it. The checks to verify that the input raster is GoogleMapsCompatible are missing and I still have to implement the test to cover the code.

CartoDB / raster-loader

Use QuadBin #45

@francois-baptiste Nice...XYZ tiles would also involve a resampling method. Would you imagine a table like: | X | Y | Z | landcat_uint8 |