developmentseed / cogeo-mosaic

Create and use COG mosaic based on mosaicJSON
https://developmentseed.org/cogeo-mosaic/
MIT License
98 stars 25 forks source link

Question: Diagnosing empty mosaic on cli create command #186

Closed zacharyDez closed 1 year ago

zacharyDez commented 2 years ago

Hi @vincentsarago,

I've been trying to test out the creation of a mosaic from a repo of radarsat imagery of the canadian federal government. I saw some issues related to s3 paths and replacing with the equivalent HTTP URLs. I can curl the URLs defined within the .txt file.

I initially created the list of images with:

aws s3 ls s3://radarsat-r1-l1-cog/1996/ --recursive | awk '{print $4}' > 1996.txt
sed -i -e 's#^#http://radarsat-r1-l1-cog.s3.amazonaws.com/#' 1996.txt

I lazily just deleted all the radarsat imagery that did not have the same configurations to avoid the datatype errors. I then ran cogeo-mosaic create 1996.txt -o mosaic.json. I'm getting the following warnings but I am not sure how to interpret them:

/home/zacharydeziel/.local/lib/python3.9/site-packages/rio_tiler/io/cogeo.py:185: UserWarning: Cannot dertermine min/max zoom based on dataset information, will default to TMS min/max zoom.
  warnings.warn(
/home/zacharydeziel/.local/lib/python3.9/site-packages/rio_tiler/io/base.py:66: UserWarning: Cannot dertermine bounds in geographic CRS, will default to (-180.0, -90.0, 180.0, 90.0).
  warnings.warn(
Get footprints  [####################################]  100%          
Get quadkey list for zoom: 0
/home/zacharydeziel/.local/lib/python3.9/site-packages/mercantile/__init__.py:77: FutureWarning: Mercantile 2.0 will require tile x and y to be within the range (0, 2 ** zoom)
  warnings.warn(
Feed Quadkey index
Iterate over quadkeys  [####################################]  100%

The resulting mosaic has no tiles:

File: mosaic.json
Backend: File
File Size: 169
Compressed: False

Profile
    MosaicJSON:       0.0.2
    Version:          1.0.0
    Name:             None
    Description:      None
    Attribution:      None

Geo
    TileMatrixSet:    WebMercatorQuad
    BoundingBox:      (-180.0, -90.0, 180.0, 90.0)
    Center:           (0.0, 0.0, 0)
    Min Zoom:         0
    Max Zoom:         24
    QuadKey Zoom:     0

Tiles
    Nb Tiles:         None
    Min Files:        None
    Max Files:        None
    Mean Files:       None

The output .txt is:

http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_C0004769_S3_19961015_182254_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_M0107472_S7_19961004_105001_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_M0202287_S7_19961018_103938_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_M0202288_S7_19961018_103952_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_M0202289_S7_19961018_104007_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_P0385952_S7_19961007_222206_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/10/RS1_P0385953_S7_19961007_222217_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_C0012334_S7_19961114_224922_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_C0014519_F5_19961120_111805_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_M0093351_EH6_19961116_113339_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_M0109175_S7_19961112_133510_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_P0393263_F1_19961115_162740_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_P0393264_F1_19961115_162746_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/11/RS1_P0451672_S2_19961103_082333_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_C0007151_S7_19961222_223828_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_C0007152_S7_19961223_112328_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_C0014851_F2_19961217_095803_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_C0018010_S3_19961213_230828_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_C0018011_S3_19961213_230836_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_M0161542_S7_19961220_220128_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_M0408657_F3_19961230_034032_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/12/RS1_P0393266_F1_19961216_162356_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/3/RS1_1024097_F5_19960305_093113_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_M0054415_S7_19960610_052834_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_M0055561_F1_19960612_223420_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_M0056541_W2_19960616_105648_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_M0106144_S7_19960612_161128_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_N0517820_S7_19960604_180642_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_P0449705_F1_19960612_223427_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_P0589220_W2_19960613_001522_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/6/RS1_P0589223_S7_19960616_002752_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_C0003543_S6_19960818_103314_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_M0060521_S4_19960801_231513_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_M0064563_S1_19960804_213347_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_M0620463_S6_19960813_142450_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_M0700301_S5_19960814_134810_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_N0408656_F4_19960811_162755_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_P0449707_F4_19960809_224250_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_P0487455_S1_19960814_114004_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_X0000001_S6_19960831_223456_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_X0000002_S7_19960815_111223_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_X0000003_F4_19960829_110520_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_X0000004_F4_19960807_223511_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/8/RS1_X0000005_F3_19960805_110526_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_C0003779_S6_19960912_214342_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_C0004698_S7_19960917_140746_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_C0004699_S7_19960926_094033_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_C0004700_S7_19960913_110303_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_M0107470_S7_19960916_125609_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_M0107473_S7_19960918_134015_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_M0109176_S7_19960922_113851_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_M0109177_S7_19960912_131113_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_M0620464_S7_19960920_141645_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_P0374972_S7_19960911_101813_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_P0374973_S7_19960911_101830_HH_SGF.tif
http://radarsat-r1-l1-cog.s3.amazonaws.com/1996/9/RS1_P0393262_F2_19960901_034049_HH_SGF.tif

Any suggestions on next steps to solve the issue? Any other tips are more than welcome.

Thanks

vincentsarago commented 2 years ago

@zacharyDez 👋

your COGs don't have CRS or geo_transform set, but use GCPS (ground control points) for geo-reference.

Sadly you can't use cogeo-mosaic cli directly for this.

  1. write a small CLI
    
    # cog_gcp_info.py
    from rio_tiler.io import GCPCOGReader

import click from rasterio.rio import options import json

@click.command() @options.file_in_arg def main(input): """Read tile.""" with GCPCOGReader(input) as cog: bounds = cog.geographic_bounds click.echo( json.dumps( { "geometry": { "type": "Polygon", "coordinates": [ [ [bounds[0], bounds[3]], [bounds[0], bounds[1]], [bounds[2], bounds[1]], [bounds[2], bounds[3]], [bounds[0], bounds[3]], ] ], }, "properties": { "path": input, "bounds": cog.geographic_bounds, "minzoom": cog.minzoom, "maxzoom": cog.maxzoom, "datatype": cog.dataset.meta["dtype"], }, "type": "Feature", } ) )

if name == 'main': main()


2. Run CLI on your files

```bash
# using parallel is optional
$ cat 1996.txt | parallel  -j 4 python -m cog_gcp_info {} > 1996.json
  1. run cogeo-mosaic cli create-from-features
$ cat 1996.json|  cogeo-mosaic create-from-features --minzoom 8 --maxzoom 12 -o 1996_mosaic.json --prop
erty path
vincentsarago commented 2 years ago

Side notes: your 255 COG, cover different zoom levels

➜  Dev cat 1996.json| jq '.properties.minzoom' | sort | uniq
10
2
3
4
5
6
7
8
9
➜  Dev cat 1996.json| jq '.properties.maxzoom' | sort | uniq
10
11
12
13
14
15
8
9

which might in fact correspond to different Product, it might be smart to first select only some scenes instead of creating one mosaic of all scenes.

Also, note that there are a lot of overlaping scenes which might impact the performance.

list of features: https://gist.github.com/vincentsarago/f099516b216b6fd7ff4ca7ac2c43f7d4 mosaic: https://gist.github.com/vincentsarago/e8896078577c9ffe9300cae7a54d31d3

vincentsarago commented 2 years ago

side note 2: the mosaic won't work natively in TiTiler, you'll need to use MosaicTilerFactory and set the dataset_reader to rio_tiler.io.GCPCOGReader

https://github.com/developmentseed/titiler/blob/master/src/titiler/mosaic/titiler/mosaic/factory.py#L44

zacharyDez commented 2 years ago

Thanks @vincentsarago ! on point as always. Will try to work through this shortly. I'll get back to you.

vincentsarago commented 2 years ago

@zacharyDez are we good to close here?

zacharyDez commented 2 years ago

yes, sorry @vincentsarago, did not get a chance to get back to this.