add_layer doesn't work with GeoTIFF file converted from xarray Dataset

[zxdawn]

I'm trying to add the layer using GeoTIFF file converted from xarray Dataset. However, it isn't shown in the map.

import numpy as np
import xarray as xr
import rioxarray
from localtileserver import TileClient, get_leaflet_tile_layer, examples
from ipyleaflet import Map

data = np.random.randint(100, size = (50, 50))
grid_coords = {'lat': np.arange(0, 50, 1), 'lon': np.arange(0, 50, 1)}

da_1 = xr.DataArray(data, dims=['lat', 'lon'], coords=grid_coords).rename('1')
da_2 = xr.DataArray(data*2, dims=['lat', 'lon'], coords=grid_coords).rename('2')

output = xr.merge([da_1, da_2])
output.rio.set_spatial_dims(x_dim='lon', y_dim='lat').rio.to_raster('test.tif')

client = TileClient('test.tif')

# Create ipyleaflet TileLayer from that server
t = get_leaflet_tile_layer(client, band=1)

# Create ipyleaflet map, add tile layer, and display
m = Map(center=client.center(), zoom=client.default_zoom)
m.add_layer(t)
m

[banesullivan]

I'm able to reproduce this... Something weird is happening here indeed. I will get to the bottom of this.

[zxdawn]

Thanks! BTW, how did you usually create GeoTIFF files for localtileserver? Maybe conversion leads to something wrong.

[giswqs]

Try leafmap numpy_to_cog. https://leafmap.org/notebooks/47_numpy_to_cog/

[banesullivan]

There's definitely a bug specifically in localtileservers' handling of this data type that needs to be addressed, but maybe using numpy_to_cog or using rio-cogeo to convert the saved TIF to a cog could mitigate this in the meantime

how did you usually create GeoTIFF files for localtileserver?

I usually use rasterio directly

[giswqs]

You can also use leafmap.image_to_cog() to convert an image to COG.
https://leafmap.org/common/#leafmap.common.image_to_cog

[zxdawn]

Thank you, guys. I tested numpy_to_cog but it still failed.

import numpy as np
import xarray as xr
from localtileserver import TileClient, get_leaflet_tile_layer
from ipyleaflet import Map
import leafmap


data = np.random.randint(100, size = (50, 50))
grid_coords = {'lat': np.arange(0, 50, 1), 'lon': np.arange(0, 50, 1)}

da_1 = xr.DataArray(data, dims=['lat', 'lon'], coords=grid_coords).rename('1')
da_2 = xr.DataArray(data*2, dims=['lat', 'lon'], coords=grid_coords).rename('2')

output = xr.merge([da_1, da_2])
leafmap.numpy_to_cog(output.to_array().values, 'test.tif')

client = TileClient('test.tif')

# Create ipyleaflet TileLayer from that server
t = get_leaflet_tile_layer(client, band=1)

# Create ipyleaflet map, add tile layer, and display
# m = Map(center=client.center(), zoom=client.default_zoom)
m = Map(center=client.center(), zoom=client.default_zoom)
m.add_layer(t)
m

[banesullivan]

I should have some free time soon to dive into this and resolve it. I have some suspicions of what's happening... stay tuned.

[zxdawn]

I realize I didn't pass the bound to numpy_to_cog . Here's the update code:

import numpy as np
import xarray as xr
from localtileserver import TileClient, get_leaflet_tile_layer
from ipyleaflet import Map
import leafmap

data = np.random.randint(100, size = (50, 50))
lat = np.arange(0, 50, 1)
lon = np.arange(0, 50, 1)
grid_coords = {'lat': lat, 'lon': lon}

da_1 = xr.DataArray(data, dims=['lat', 'lon'], coords=grid_coords).rename('1')
da_2 = xr.DataArray(data*2, dims=['lat', 'lon'], coords=grid_coords).rename('2')

output = xr.merge([da_1, da_2])

leafmap.numpy_to_cog(output.to_array().values, 'test.tif', bounds=(lon.min(), lat.min(), lon.max(),lat.max()), dtype='float32')

client = TileClient('test.tif')

# Create ipyleaflet TileLayer from that server
t = get_leaflet_tile_layer(client, band=1)

# Create ipyleaflet map, add tile layer, and display
m = Map(center=client.center(), zoom=2)#client.default_zoom)
m.add_layer(t)
m

I can see the data shown now. However, it only works with zoom=2. If I zoom in or out, the tile will disappear.

[banesullivan]

The issue with it only appearing at certain zoom levels is separate from the original issue. To fix this, you need to increase both max_native_zoom and max_zoom when using get_leaflet_tile_layer(). I actually had this same issue the other day and pushed a fix that isn't release in ff70753

As a quick fix, use max_zoom=30 and max_native_zoom=30

[giswqs]

In leafmap, I already fixed this. You should be able to use it directly without setting these two parameters.
https://github.com/opengeos/leafmap/blob/master/leafmap/common.py#L2813

[banesullivan]

This will be fixed in the next release (0.10.0) given the switch to the more robust rio-tiler backend

(verified example in #159 (comment))