geopandas / contextily

Context geo-tiles in Python
https://contextily.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
530 stars 83 forks source link

Reset_extent is not behaving as expected with local file when raster extents is smaller than the shape extents. #213

Open cordmaur opened 1 year ago

cordmaur commented 1 year ago

Hi,

I've made some tests to illustrate the problem I am facing. Here is my initial shape (world) and a syntetic raster with a square in Minas Gerais, Brazil. Visualizing it with reset_extent set to False: image

If we filter the shape to Brazil (smaller than the raster) and set reset_extent to True, it will visualize correctly the raster. image

But if we leave reset_extent = True and tries to visualize the world, it will stretch the raster to the world extents and it has no meaning at all... image

Am I missing something here?

Thanks, Mauricio

darribas commented 1 year ago

mmm... hard to tell without the underlying data. Is it possible at all that the crs differ between the two? My first instict would be to set crs=world.crs to be sure they're both aligned, but maybe chuva.crs is the same?

cordmaur commented 1 year ago

Hi @darribas , here is a code snippet that pulls the raster (small) from the internet. Besides the problem with the extents, I've noticed that the result with reset_extent=True is also upside-down. :)

import geopandas as gpd
import xarray as xr
import rioxarray as xrio
import contextily as cx

# get the world from geopandas
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
print(world.crs)

# get rain from INPE
# !curl -O http://ftp.cptec.inpe.br/modelos/tempo/MERGE/GPM/CLIMATOLOGY/MONTHLY_ACCUMULATED_YEARLY/MERGE_CPTEC_acum_sep_2022.nc
rain = xr.open_dataset('./MERGE_CPTEC_acum_sep_2022.nc')['pacum'].squeeze()

# set the CRS as the same from the world (EPSG:4326)
rain = rain.rio.write_crs(world.crs)

# save Rain as geotiff to be used in Contextily
rain = rain.rename({'lat': 'y', 'lon': 'x'})
rain.rio.to_raster('./rain.tif', compress='deflate')

# Open the rain as tiff just for checking
rain_tif = xrio.open_rasterio('./rain.tif')
rain_tif.plot(aspect='equal', figsize=(5, 5))

image

ax = world.plot(facecolor='none', figsize=(10, 7))
cx.add_basemap(ax=ax, source='./rain.tif', reset_extent=False)

image

rain.rio.to_raster('./rain.tif', compress='deflate')
ax = world.plot(facecolor='none', figsize=(10, 7))

cx.add_basemap(ax=ax, source='./rain.tif', reset_extent=True)

image