Spilhaus from PROJ

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.11%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22marimo%3E%3D0.24%22%2C%0A%23%20%20%20%20%20%22numpy%22%2C%0A%23%20%20%20%20%20%22pyproj%3E%3D3.7%22%2C%0A%23%20%20%20%20%20%22rasterio%3E%3D1.4%22%2C%0A%23%20%20%20%20%20%22pyogrio%3E%3D0.13%22%2C%0A%23%20%20%20%20%20%22geopandas%3E%3D1.0%22%2C%0A%23%20%20%20%20%20%22shapely%3E%3D2.0%22%2C%0A%23%20%20%20%20%20%22matplotlib%3E%3D3.8%22%2C%0A%23%20%5D%0A%23%20%2F%2F%2F%0A%22%22%22%0ASpilhaus%20from%20PROJ%2C%20through%20the%20Python%20stack.%0A%0AThe%20Spilhaus%20projection%20landed%20in%20PROJ%209.6.0%20(March%202025).%20This%20notebook%0Abuilds%20a%20world-ocean%20map%20with%20it%20using%20only%20public%20sources%20read%20remotely%0A(%2Fvsicurl)%2C%20and%20-%20the%20real%20point%20-%20reports%20which%20PROJ%20each%20downstream%0Apackage%20actually%20carries%2C%20because%20in%20a%20pip%20environment%20every%20package%0Athat%20bundles%20GDAL%20or%20PROJ%20bundles%20its%20own%20copy.%0A%0ARun%20it%20yourself%3A%0A%0A%20%20%20%20uvx%20marimo%20edit%20--sandbox%20spilhaus.py%0A%0AEnvironment%20overrides%3A%0A%0A%20%20%20%20SPILHAUS_RASTER%20%20%20any%20GDAL-readable%20global%20raster%20(default%3A%20GEBCO%202024%20COG)%0A%20%20%20%20SPILHAUS_OCEAN%20%20%20%20any%20OGR-readable%20ocean%20polygon%20layer%20(default%3A%20Natural%20Earth%20110m)%0A%20%20%20%20SPILHAUS_NPIX%20%20%20%20%20output%20grid%20size%20in%20pixels%20(default%201600)%0A%22%22%22%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.24.2%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22%2C%20app_title%3D%22Spilhaus%20from%20PROJ%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3Csmall%3EAlso%20in%20this%20series%3A%0A%20%20%20%20%3Ca%20href%3D%22..%2Fr%2F%22%3Ethe%20R%20version%3C%2Fa%3E%20(gdalraster%2C%20wk%2C%20ximage%3B%20rendered%20on%20two%20Ubuntu%20runners)%20%7C%0A%20%20%20%20%3Ca%20href%3D%22..%2Fwasm%2F%22%3Ethe%20browser%20version%3C%2Fa%3E%20(pyproj%20in%20Pyodide)%20%7C%0A%20%20%20%20%3Ca%20href%3D%22..%2F%22%3Eindex%3C%2Fa%3E%20%7C%20%3Ca%20href%3D%22https%3A%2F%2Fgithub.com%2Fmdsumner%2Fspilhaus-tracker%22%3Esource%3C%2Fa%3E%20%7C%0A%20%20%20%20%3Ca%20href%3D%22https%3A%2F%2Fwww.hypertidy.org%2Fposts%2F2026-09-20_road-to-spilhaus%2F%22%3Ethe%20story%20so%20far%3C%2Fa%3E%3C%2Fsmall%3E%0A%0A%20%20%20%20%23%20Spilhaus%2C%20as%20it%20flows%20downstream%20from%20PROJ%0A%0A%20%20%20%20%60%2Bproj%3Dspilhaus%60%20was%20added%20in%20**PROJ%209.6.0**%20(March%202025).%20Whether%20*you*%0A%20%20%20%20can%20use%20it%20from%20Python%20depends%20not%20on%20what%20PROJ%20is%20installed%20on%20your%0A%20%20%20%20machine%2C%20but%20on%20which%20PROJ%20each%20package%20carries%20with%20it.%20In%20a%20%60pip%60%0A%20%20%20%20environment%2C%20%60pyproj%60%2C%20%60rasterio%60%2C%20%60pyogrio%60%20and%20the%20%60osgeo%60%20bindings%20each%0A%20%20%20%20ship%20(or%20link)%20a%20**separate**%20copy%20of%20PROJ%20and%20GDAL.%0A%0A%20%20%20%20This%20notebook%3A%0A%0A%20%20%20%201.%20tabulates%20the%20PROJ%20and%20GDAL%20behind%20each%20package%20and%20probes%20each%20one%20for%20%60%2Bproj%3Dspilhaus%60%3B%0A%20%20%20%202.%20builds%20a%20Spilhaus%20world-ocean%20map%20from%20public%20remote%20sources%20(bathymetry%20COG%20over%20%60%2Fvsicurl%60%2C%0A%20%20%20%20%20%20%20Natural%20Earth%20ocean%20polygons%2C%20a%20handful%20of%20hand-typed%20landmarks)%3B%0A%20%20%20%203.%20uses%20whichever%20library%20*can*%20do%20the%20transform%2C%20and%20says%20which%20one%20it%20was.%0A%0A%20%20%20%20Nothing%20is%20downloaded%20up%20front%3B%20every%20source%20is%20read%20remotely%20via%20GDAL's%20virtual%20file%20systems.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20os%0A%20%20%20%20import%20shutil%0A%20%20%20%20import%20subprocess%0A%20%20%20%20import%20sys%0A%20%20%20%20import%20warnings%0A%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20warnings.filterwarnings(%22ignore%22%2C%20category%3DFutureWarning)%0A%0A%20%20%20%20%23%20Remote-read%20hygiene%20for%20%2Fvsicurl%3A%20no%20directory%20listings%2C%20only%20fetch%20what%20is%20asked%20for.%0A%20%20%20%20os.environ.setdefault(%22GDAL_DISABLE_READDIR_ON_OPEN%22%2C%20%22EMPTY_DIR%22)%0A%20%20%20%20os.environ.setdefault(%22GDAL_HTTP_MULTIRANGE%22%2C%20%22YES%22)%0A%20%20%20%20os.environ.setdefault(%22GDAL_HTTP_MERGE_CONSECUTIVE_RANGES%22%2C%20%22YES%22)%0A%20%20%20%20os.environ.setdefault(%22VSI_CACHE%22%2C%20%22TRUE%22)%0A%20%20%20%20os.environ.setdefault(%22VSI_CACHE_SIZE%22%2C%20str(64%20*%201024%20*%201024))%0A%0A%20%20%20%20SPILHAUS%20%3D%20%22%2Bproj%3Dspilhaus%22%0A%0A%20%20%20%20%23%20The%20same%20map%20with%20every%20parameter%20spelled%20out%20(values%20from%20the%20PROJ%20docs).%0A%20%20%20%20%23%20Bare%20%22%2Bproj%3Dspilhaus%22%20relies%20on%20defaults%20that%20live%20in%20PROJ's%20source%20code.%0A%20%20%20%20SPILHAUS_EXPLICIT%20%3D%20(%0A%20%20%20%20%20%20%20%20%22%2Bproj%3Dspilhaus%20%2Blon_0%3D66.94970198%20%2Blat_0%3D-49.56371678%20%22%0A%20%20%20%20%20%20%20%20%22%2Bazi%3D40.17823482%20%2Brot%3D45%20%2Bk_0%3D1%20%2Bx_0%3D0%20%2By_0%3D0%20%2BR%3D6378137%22%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20A%20different%20member%20of%20the%20same%20family%3A%20naive%20recentring%20on%20the%20Atlantic.%0A%20%20%20%20%23%20Same%20projection%20method%2C%20different%20parameters%2C%20and%20the%20cut%20now%20runs%20through%20ocean.%0A%20%20%20%20SPILHAUS_ALT%20%3D%20%22%2Bproj%3Dspilhaus%20%2Blon_0%3D-30%20%2Blat_0%3D-20%20%2Bazi%3D20%20%2Brot%3D45%22%0A%0A%20%20%20%20GEBCO_URL%20%3D%20%22%2Fvsicurl%2Fhttps%3A%2F%2Fdata.source.coop%2Falexgleith%2Fgebco-2024%2FGEBCO_2024.tif%22%0A%20%20%20%20RASTER_URL%20%3D%20os.environ.get(%22SPILHAUS_RASTER%22%2C%20GEBCO_URL)%0A%0A%20%20%20%20NE_OCEAN_URL%20%3D%20(%0A%20%20%20%20%20%20%20%20%22%2Fvsicurl%2Fhttps%3A%2F%2Fraw.githubusercontent.com%2Fnvkelso%2Fnatural-earth-vector%2F%22%0A%20%20%20%20%20%20%20%20%22master%2F110m_physical%2Fne_110m_ocean.shp%22%0A%20%20%20%20)%0A%20%20%20%20OCEAN_URL%20%3D%20os.environ.get(%22SPILHAUS_OCEAN%22%2C%20NE_OCEAN_URL)%0A%0A%20%20%20%20NPIX%20%3D%20int(os.environ.get(%22SPILHAUS_NPIX%22%2C%20%221600%22))%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20NPIX%2C%0A%20%20%20%20%20%20%20%20OCEAN_URL%2C%0A%20%20%20%20%20%20%20%20RASTER_URL%2C%0A%20%20%20%20%20%20%20%20SPILHAUS%2C%0A%20%20%20%20%20%20%20%20SPILHAUS_ALT%2C%0A%20%20%20%20%20%20%20%20SPILHAUS_EXPLICIT%2C%0A%20%20%20%20%20%20%20%20np%2C%0A%20%20%20%20%20%20%20%20os%2C%0A%20%20%20%20%20%20%20%20shutil%2C%0A%20%20%20%20%20%20%20%20subprocess%2C%0A%20%20%20%20%20%20%20%20sys%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%201.%20Who%20has%20which%20PROJ%3F%0A%0A%20%20%20%20Each%20row%20is%20one%20package's%20view%20of%20the%20world.%20The%20final%20column%20is%20the%0A%20%20%20%20result%20of%20asking%20that%20package%20to%20create%20a%20CRS%20from%20%60%2Bproj%3Dspilhaus%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(SPILHAUS%2C%20mo%2C%20shutil%2C%20subprocess%2C%20sys)%3A%0A%20%20%20%20import%20importlib.metadata%20as%20_md%0A%0A%20%20%20%20def%20_ver(pkg)%3A%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20_md.version(pkg)%0A%20%20%20%20%20%20%20%20except%20_md.PackageNotFoundError%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20None%0A%0A%20%20%20%20probes%20%3D%20%5B%5D%20%20%23%20rows%20for%20the%20table%0A%20%20%20%20engines%20%3D%20%7B%7D%20%20%23%20name%20-%3E%20factory(proj_string)%20-%3E%20dict(forward%3Dcallable%2C%20inverse%3Dcallable)%0A%0A%20%20%20%20%23%20---%20pyproj%20-------------------------------------------------------------%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20import%20pyproj%0A%0A%20%20%20%20%20%20%20%20_row%20%3D%20dict(package%3D%22pyproj%22%2C%20version%3Dpyproj.__version__%2C%20gdal%3D%22-%22%2C%20proj%3Dpyproj.proj_version_str)%0A%20%20%20%20%20%20%20%20def%20_pyproj_factory(ps%2C%20_pp%3Dpyproj)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20crs%20%3D%20_pp.CRS.from_proj4(ps)%0A%20%20%20%20%20%20%20%20%20%20%20%20fwd%20%3D%20_pp.Transformer.from_crs(4326%2C%20crs%2C%20always_xy%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20inv%20%3D%20_pp.Transformer.from_crs(crs%2C%204326%2C%20always_xy%3DTrue)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20dict(forward%3Dfwd.transform%2C%20inverse%3Dinv.transform)%0A%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_pyproj_factory(SPILHAUS)%20%20%23%20the%20probe%0A%20%20%20%20%20%20%20%20%20%20%20%20engines%5B%22pyproj%22%5D%20%3D%20_pyproj_factory%0A%20%20%20%20%20%20%20%20%20%20%20%20_row%5B%22spilhaus%22%5D%20%3D%20%22ok%22%0A%20%20%20%20%20%20%20%20except%20Exception%20as%20e%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20_row%5B%22spilhaus%22%5D%20%3D%20f%22FAIL%3A%20%7Btype(e).__name__%7D%22%0A%20%20%20%20%20%20%20%20probes.append(_row)%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20probes.append(dict(package%3D%22pyproj%22%2C%20version%3DNone%2C%20gdal%3DNone%2C%20proj%3DNone%2C%20spilhaus%3D%22not%20installed%22))%0A%0A%20%20%20%20%23%20---%20rasterio%20-----------------------------------------------------------%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20import%20rasterio%0A%20%20%20%20%20%20%20%20import%20rasterio.crs%0A%20%20%20%20%20%20%20%20from%20rasterio.warp%20import%20transform%20as%20_rio_transform%0A%0A%20%20%20%20%20%20%20%20_proj%20%3D%20getattr(rasterio%2C%20%22__proj_version__%22%2C%20None)%20or%20%22%3F%22%0A%20%20%20%20%20%20%20%20_row%20%3D%20dict(package%3D%22rasterio%22%2C%20version%3Drasterio.__version__%2C%20gdal%3Drasterio.__gdal_version__%2C%20proj%3Dstr(_proj))%0A%20%20%20%20%20%20%20%20def%20_rasterio_factory(ps%2C%20_rio%3Drasterio%2C%20_t%3D_rio_transform)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20crs%20%3D%20_rio.crs.CRS.from_proj4(ps)%0A%20%20%20%20%20%20%20%20%20%20%20%20wgs%20%3D%20_rio.crs.CRS.from_epsg(4326)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20dict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20forward%3Dlambda%20x%2C%20y%3A%20_t(wgs%2C%20crs%2C%20list(x)%2C%20list(y))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20inverse%3Dlambda%20x%2C%20y%3A%20_t(crs%2C%20wgs%2C%20list(x)%2C%20list(y))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_rasterio_factory(SPILHAUS)%20%20%23%20the%20probe%0A%20%20%20%20%20%20%20%20%20%20%20%20engines%5B%22rasterio%22%5D%20%3D%20_rasterio_factory%0A%20%20%20%20%20%20%20%20%20%20%20%20_row%5B%22spilhaus%22%5D%20%3D%20%22ok%22%0A%20%20%20%20%20%20%20%20except%20Exception%20as%20e%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20_row%5B%22spilhaus%22%5D%20%3D%20f%22FAIL%3A%20%7Btype(e).__name__%7D%22%0A%20%20%20%20%20%20%20%20probes.append(_row)%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20probes.append(dict(package%3D%22rasterio%22%2C%20version%3DNone%2C%20gdal%3DNone%2C%20proj%3DNone%2C%20spilhaus%3D%22not%20installed%22))%0A%0A%20%20%20%20%23%20---%20pyogrio%20(vector%20I%2FO%3B%20links%20GDAL%2C%20exposes%20no%20CRS-transform%20API)%20------%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20import%20pyogrio%0A%0A%20%20%20%20%20%20%20%20probes.append(%0A%20%20%20%20%20%20%20%20%20%20%20%20dict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20package%3D%22pyogrio%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20version%3Dpyogrio.__version__%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20gdal%3Dpyogrio.__gdal_version_string__%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20proj%3D%22(static%20in%20libgdal)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20spilhaus%3D%22n%2Fa%20(I%2FO%20only)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20probes.append(dict(package%3D%22pyogrio%22%2C%20version%3DNone%2C%20gdal%3DNone%2C%20proj%3DNone%2C%20spilhaus%3D%22not%20installed%22))%0A%0A%20%20%20%20%23%20---%20osgeo%20(GDAL's%20own%20Python%20bindings)%20---------------------------------%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20from%20osgeo%20import%20gdal%20as%20_gdal%0A%20%20%20%20%20%20%20%20from%20osgeo%20import%20osr%20as%20_osr%0A%0A%20%20%20%20%20%20%20%20_osr.UseExceptions()%0A%20%20%20%20%20%20%20%20_gdal.UseExceptions()%0A%20%20%20%20%20%20%20%20_proj%20%3D%20f%22%7B_osr.GetPROJVersionMajor()%7D.%7B_osr.GetPROJVersionMinor()%7D.%7B_osr.GetPROJVersionMicro()%7D%22%0A%20%20%20%20%20%20%20%20_row%20%3D%20dict(package%3D%22osgeo%20(GDAL%20bindings)%22%2C%20version%3D_gdal.__version__%2C%20gdal%3D_gdal.__version__%2C%20proj%3D_proj)%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_srs%20%3D%20_osr.SpatialReference()%0A%20%20%20%20%20%20%20%20%20%20%20%20_srs.ImportFromProj4(SPILHAUS)%0A%20%20%20%20%20%20%20%20%20%20%20%20_row%5B%22spilhaus%22%5D%20%3D%20%22ok%22%0A%20%20%20%20%20%20%20%20except%20Exception%20as%20e%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20_row%5B%22spilhaus%22%5D%20%3D%20f%22FAIL%3A%20%7Btype(e).__name__%7D%22%0A%20%20%20%20%20%20%20%20probes.append(_row)%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20probes.append(dict(package%3D%22osgeo%20(GDAL%20bindings)%22%2C%20version%3DNone%2C%20gdal%3DNone%2C%20proj%3DNone%2C%20spilhaus%3D%22not%20installed%22))%0A%0A%20%20%20%20%23%20---%20geopandas%20%2F%20shapely%3A%20ride%20on%20pyproj%20and%20GEOS%20------------------------%0A%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20import%20geopandas%20as%20_gpd%0A%20%20%20%20%20%20%20%20import%20shapely%20as%20_shp%0A%0A%20%20%20%20%20%20%20%20probes.append(%0A%20%20%20%20%20%20%20%20%20%20%20%20dict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20package%3D%22geopandas%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20version%3D_gpd.__version__%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20gdal%3D%22(via%20pyogrio)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20proj%3D%22(via%20pyproj)%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20spilhaus%3D%22inherits%20pyproj%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20probes.append(%0A%20%20%20%20%20%20%20%20%20%20%20%20dict(package%3D%22shapely%22%2C%20version%3D_shp.__version__%2C%20gdal%3D%22-%22%2C%20proj%3Df%22GEOS%20%7B_shp.geos_version_string%7D%22%2C%20spilhaus%3D%22n%2Fa%22)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20except%20ImportError%3A%0A%20%20%20%20%20%20%20%20pass%0A%0A%20%20%20%20%23%20---%20command%20line%20tools%20on%20PATH%20-------------------------------------------%0A%20%20%20%20def%20_cli_version(cmd%2C%20args)%3A%0A%20%20%20%20%20%20%20%20exe%20%3D%20shutil.which(cmd)%0A%20%20%20%20%20%20%20%20if%20not%20exe%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20None%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20out%20%3D%20subprocess.run(%5Bexe%2C%20*args%5D%2C%20capture_output%3DTrue%2C%20text%3DTrue%2C%20timeout%3D30)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20(out.stdout%20or%20out.stderr).strip().splitlines()%5B0%5D%0A%20%20%20%20%20%20%20%20except%20Exception%20as%20e%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20f%22error%3A%20%7Be%7D%22%0A%0A%20%20%20%20cli%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22gdalinfo%20--version%22%3A%20_cli_version(%22gdalinfo%22%2C%20%5B%22--version%22%5D)%2C%0A%20%20%20%20%20%20%20%20%22gdal%20--version%20(unified%20CLI%2C%20GDAL%20%3E%3D%203.11)%22%3A%20_cli_version(%22gdal%22%2C%20%5B%22--version%22%5D)%2C%0A%20%20%20%20%20%20%20%20%22projinfo%20--version%22%3A%20_cli_version(%22projinfo%22%2C%20%5B%22--version%22%5D)%2C%0A%20%20%20%20%20%20%20%20%22rio%20--version%22%3A%20_cli_version(%22rio%22%2C%20%5B%22--version%22%5D)%2C%0A%20%20%20%20%7D%0A%0A%20%20%20%20capable%20%3D%20%5Bk%20for%20k%20in%20engines%5D%0A%20%20%20%20mo.vstack(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.md(f%22Python%20%60%7Bsys.version.split()%5B0%5D%7D%60%20at%20%60%7Bsys.executable%7D%60%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.ui.table(probes%2C%20selection%3DNone%2C%20label%3D%22Native%20library%20provenance%2C%20one%20row%20per%20package%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22**Command-line%20tools%20on%20PATH%3A**%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%2B%20%22%5Cn%22.join(f%22-%20%60%7Bk%7D%60%3A%20%7Bv%20if%20v%20else%20'not%20found'%7D%22%20for%20k%2C%20v%20in%20cli.items())%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22Libraries%20that%20can%20build%20%60%2Bproj%3Dspilhaus%60%20here%3A%20**%7B'%2C%20'.join(capable)%20or%20'none'%7D**.%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22The%20rest%20of%20this%20notebook%20uses%20the%20first%20of%20these%20for%20coordinate%20transforms.%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20kind%3D%22success%22%20if%20capable%20else%20%22danger%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20return%20capable%2C%20engines%2C%20probes%0A%0A%0A%40app.cell%0Adef%20_(SPILHAUS%2C%20capable%2C%20engines%2C%20np)%3A%0A%20%20%20%20%23%20Forward%2Finverse%20constructors%20for%20any%20proj%20string%2C%20from%20whichever%20library%20was%20capable.%0A%20%20%20%20if%20not%20capable%3A%0A%20%20%20%20%20%20%20%20raise%20RuntimeError(%22No%20installed%20library%20can%20create%20%2Bproj%3Dspilhaus%3B%20nothing%20below%20can%20run.%22)%0A%0A%20%20%20%20ENGINE%20%3D%20capable%5B0%5D%0A%0A%20%20%20%20def%20make_forward(proj_string)%3A%0A%20%20%20%20%20%20%20%20_fwd%20%3D%20engines%5BENGINE%5D(proj_string)%5B%22forward%22%5D%0A%0A%20%20%20%20%20%20%20%20def%20forward(lon%2C%20lat)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22lon%2Flat%20arrays%20-%3E%20projected%20x%2Fy%20arrays%20(NaN%20where%20undefined).%22%22%22%0A%20%20%20%20%20%20%20%20%20%20%20%20lon%20%3D%20np.asarray(lon%2C%20dtype%3Dfloat).ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20lat%20%3D%20np.asarray(lat%2C%20dtype%3Dfloat).ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20x%2C%20y%20%3D%20_fwd(lon.tolist()%2C%20lat.tolist())%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20np.asarray(x%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20%20%20%20%20y%20%3D%20np.asarray(y%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20%20%20%20%20bad%20%3D%20~np.isfinite(x)%20%7C%20~np.isfinite(y)%20%7C%20(np.abs(x)%20%3E%201e15)%20%7C%20(np.abs(y)%20%3E%201e15)%0A%20%20%20%20%20%20%20%20%20%20%20%20x%5Bbad%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20y%5Bbad%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20x%2C%20y%0A%0A%20%20%20%20%20%20%20%20return%20forward%0A%0A%20%20%20%20def%20make_inverse(proj_string)%3A%0A%20%20%20%20%20%20%20%20_inv%20%3D%20engines%5BENGINE%5D(proj_string)%5B%22inverse%22%5D%0A%0A%20%20%20%20%20%20%20%20def%20inverse(x%2C%20y)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%22%22projected%20x%2Fy%20-%3E%20lon%2Flat%2C%20one%20point%20at%20a%20time%20so%20singular%20corners%20fail%20individually.%22%22%22%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20np.asarray(x%2C%20dtype%3Dfloat).ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20y%20%3D%20np.asarray(y%2C%20dtype%3Dfloat).ravel()%0A%20%20%20%20%20%20%20%20%20%20%20%20lon%20%3D%20np.full_like(x%2C%20np.nan)%0A%20%20%20%20%20%20%20%20%20%20%20%20lat%20%3D%20np.full_like(y%2C%20np.nan)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20range(len(x))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20a%2C%20b%20%3D%20_inv(%5Bx%5Bi%5D%5D%2C%20%5By%5Bi%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20lon%5Bi%5D%2C%20lat%5Bi%5D%20%3D%20float(a%5B0%5D)%2C%20float(b%5B0%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20except%20Exception%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20pass%0A%20%20%20%20%20%20%20%20%20%20%20%20bad%20%3D%20~np.isfinite(lon)%20%7C%20(np.abs(lon)%20%3E%20360)%20%7C%20(np.abs(lat)%20%3E%2090)%0A%20%20%20%20%20%20%20%20%20%20%20%20lon%5Bbad%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20lat%5Bbad%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20lon%2C%20lat%0A%0A%20%20%20%20%20%20%20%20return%20inverse%0A%0A%20%20%20%20forward%20%3D%20make_forward(SPILHAUS)%0A%20%20%20%20inverse%20%3D%20make_inverse(SPILHAUS)%0A%20%20%20%20return%20ENGINE%2C%20forward%2C%20inverse%2C%20make_forward%0A%0A%0A%40app.cell%0Adef%20_(ENGINE%2C%20mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20rf%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%202.%20The%20projection's%20domain%0A%0A%20%20%20%20%20%20%20%20Spilhaus%20maps%20the%20sphere%20onto%20a%20square.%20PROJ%20does%20not%20publish%20the%20square's%0A%20%20%20%20%20%20%20%20half-width%20as%20metadata%2C%20so%20we%20forward-project%20a%20dense%20lon%2Flat%20grid%20(with%0A%20%20%20%20%20%20%20%20**%7BENGINE%7D**)%20and%20take%20the%20extent.%20The%20square's%20boundary%20is%20the%20*cut*%3A%20the%0A%20%20%20%20%20%20%20%20line%20along%20which%20the%20sphere%20is%20opened%2C%20chosen%20to%20run%20almost%20entirely%20through%0A%20%20%20%20%20%20%20%20land%20so%20the%20ocean%20is%20one%20connected%20piece.%20Inverse-projecting%20the%20boundary%0A%20%20%20%20%20%20%20%20back%20to%20lon%2Flat%20traces%20the%20cut%20on%20the%20globe.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(forward%2C%20inverse%2C%20np)%3A%0A%20%20%20%20def%20square_half(fwd%2C%20n%3D1441)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Half-width%20of%20the%20projected%20square%3A%20forward-project%20a%20dense%20grid%20and%20take%20the%20extent.%22%22%22%0A%20%20%20%20%20%20%20%20lon%2C%20lat%20%3D%20np.meshgrid(np.linspace(-180%2C%20180%2C%20n)%2C%20np.linspace(-90%2C%2090%2C%20(n%20%2B%201)%20%2F%2F%202))%0A%20%20%20%20%20%20%20%20x%2C%20y%20%3D%20fwd(lon%2C%20lat)%0A%20%20%20%20%20%20%20%20half%20%3D%20float(np.nanmax(np.abs(np.concatenate(%5Bx%2C%20y%5D))))%0A%20%20%20%20%20%20%20%20%23%20the%20true%20half-width%20is%20slightly%20larger%20than%20any%20sampled%20point%0A%20%20%20%20%20%20%20%20return%20float(np.ceil(half%20%2F%201000.0)%20*%201000.0)%0A%0A%20%20%20%20HALF%20%3D%20square_half(forward)%0A%0A%20%20%20%20%23%20Trace%20the%20cut%3A%20inverse-project%20the%20square%20boundary%20just%20inside%20the%20edge.%0A%20%20%20%20_n%20%3D%20300%0A%20%20%20%20_t%20%3D%20np.linspace(-1%2C%201%2C%20_n)%0A%20%20%20%20_edge%20%3D%200.996%20*%20HALF%0A%20%20%20%20_boundary%20%3D%20np.concatenate(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20np.c_%5B_t%20*%20_edge%2C%20np.full(_n%2C%20-_edge)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.c_%5Bnp.full(_n%2C%20_edge)%2C%20_t%20*%20_edge%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.c_%5B_t%5B%3A%3A-1%5D%20*%20_edge%2C%20np.full(_n%2C%20_edge)%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20np.c_%5Bnp.full(_n%2C%20-_edge)%2C%20_t%5B%3A%3A-1%5D%20*%20_edge%5D%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20cut_lon%2C%20cut_lat%20%3D%20inverse(_boundary%5B%3A%2C%200%5D%2C%20_boundary%5B%3A%2C%201%5D)%0A%0A%20%20%20%20%23%20the%20two%20singular%20corner%20points%2C%20from%20the%20literature%3A%20(115E%2C%2030N)%20and%20(65W%2C%2030S)%0A%20%20%20%20corners%20%3D%20%7B%22corner%20A%20(Asia)%22%3A%20(115.0%2C%2030.0)%2C%20%22corner%20B%20(South%20America)%22%3A%20(-65.0%2C%20-30.0)%7D%0A%20%20%20%20return%20HALF%2C%20corners%2C%20cut_lat%2C%20cut_lon%2C%20square_half%0A%0A%0A%40app.cell%0Adef%20_(HALF%2C%20corners%2C%20cut_lat%2C%20cut_lon%2C%20mo%2C%20np)%3A%0A%20%20%20%20import%20matplotlib%0A%0A%20%20%20%20matplotlib.use(%22Agg%22)%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%204.2))%0A%20%20%20%20_ax.scatter(cut_lon%2C%20cut_lat%2C%20s%3D4%2C%20color%3D%22%23b8342c%22%2C%20label%3D%22square%20boundary%2C%20inverse-projected%20(the%20cut)%22)%0A%20%20%20%20for%20_name%2C%20(_lo%2C%20_la)%20in%20corners.items()%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_lo%2C%20_la%2C%20%22o%22%2C%20ms%3D7%2C%20mfc%3D%22none%22%2C%20mec%3D%22black%22%2C%20mew%3D1.2)%0A%20%20%20%20%20%20%20%20_ax.annotate(_name%2C%20(_lo%2C%20_la)%2C%20xytext%3D(6%2C%206)%2C%20textcoords%3D%22offset%20points%22%2C%20fontsize%3D8)%0A%20%20%20%20_ax.set_xlim(-180%2C%20180)%0A%20%20%20%20_ax.set_ylim(-90%2C%2090)%0A%20%20%20%20_ax.set_aspect(%22equal%22)%0A%20%20%20%20_ax.grid(True%2C%20lw%3D0.3%2C%20alpha%3D0.6)%0A%20%20%20%20_ax.set_title(f%22Where%20the%20Spilhaus%20cut%20falls%20on%20the%20globe%20%20(square%20half-width%20%7BHALF%2F1e6%3A.3f%7D%20Mm)%22)%0A%20%20%20%20_ax.legend(loc%3D%22lower%20left%22%2C%20fontsize%3D8)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20mo.vstack(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22Half-width%20of%20the%20square%3A%20**%7BHALF%3A%2C.0f%7D%20m**%3B%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22%7Bnp.isfinite(cut_lon).sum()%7D%20of%20%7Blen(cut_lon)%7D%20boundary%20samples%20inverted%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22(the%20misses%20cluster%20at%20the%20two%20singular%20corners).%22%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20_fig%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20return%20matplotlib%2C%20plt%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%203.%20Raster%3A%20bathymetry%20over%20%60%2Fvsicurl%60%0A%0A%20%20%20%20Default%20source%20is%20Alex%20Leith's%20GEBCO%202024%20COG%20on%20source.coop%2C%20read%20with%0A%20%20%20%20rasterio's%20%60WarpedVRT%60%20-%20GDAL's%20warper%2C%20driven%20from%20Python.%20We%20open%20the%0A%20%20%20%20source%20at%20the%20overview%20level%20closest%20to%20the%20output%20resolution%20so%20the%0A%20%20%20%204%20GB%20file%20is%20never%20read%20in%20full%3B%20only%20the%20tiles%20needed%20come%20over%20the%20wire.%0A%0A%20%20%20%20%3E%20Set%20%60SPILHAUS_RASTER%60%20to%20point%20at%20another%20global%20raster%20(any%20GDAL%20path%2C%0A%20%20%20%20%3E%20%60%2Fvsicurl%2F%60%2C%20%60%2Fvsis3%2F%60%2C%20a%20WMTS%20descriptor%2C%20...).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(HALF%2C%20NPIX%2C%20RASTER_URL%2C%20SPILHAUS%2C%20mo)%3A%0A%20%20%20%20import%20time%0A%0A%20%20%20%20import%20rasterio%20as%20rio%0A%20%20%20%20from%20rasterio.enums%20import%20Resampling%0A%20%20%20%20from%20rasterio.transform%20import%20from_bounds%0A%20%20%20%20from%20rasterio.vrt%20import%20WarpedVRT%0A%0A%20%20%20%20def%20_pick_overview(src%2C%20dst_res_m)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Index%20into%20src.overviews(1)%20whose%20resolution%20is%20just%20finer%20than%20the%20target%2C%20else%20None.%22%22%22%0A%20%20%20%20%20%20%20%20ovr%20%3D%20src.overviews(1)%0A%20%20%20%20%20%20%20%20if%20not%20ovr%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20None%0A%20%20%20%20%20%20%20%20%23%20degrees%20per%20source%20pixel%20at%20level%20k%2C%20converted%20to%20metres%20at%20the%20equator%0A%20%20%20%20%20%20%20%20base_m%20%3D%20abs(src.transform.a)%20*%20111_320.0%0A%20%20%20%20%20%20%20%20levels%20%3D%20%5B(k%2C%20base_m%20*%20f)%20for%20k%2C%20f%20in%20enumerate(ovr)%5D%0A%20%20%20%20%20%20%20%20finer%20%3D%20%5Bk%20for%20k%2C%20res%20in%20levels%20if%20res%20%3C%3D%20dst_res_m%20%2F%201.5%5D%0A%20%20%20%20%20%20%20%20return%20max(finer)%20if%20finer%20else%20None%0A%0A%20%20%20%20def%20warp_square(proj_string%2C%20half%2C%20npix%2C%20url%3DRASTER_URL)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Warp%20the%20source%20raster%20onto%20the%20%5B-half%2C%20half%5D%5E2%20square%20of%20proj_string%20at%20npix%20x%20npix.%0A%0A%20%20%20%20%20%20%20%20Returns%20(array%20with%20alpha%20band%20last%2C%20metadata%20dict).%22%22%22%0A%20%20%20%20%20%20%20%20dst_crs%20%3D%20rio.crs.CRS.from_proj4(proj_string)%0A%20%20%20%20%20%20%20%20dst_transform%20%3D%20from_bounds(-half%2C%20-half%2C%20half%2C%20half%2C%20npix%2C%20npix)%0A%20%20%20%20%20%20%20%20res_m%20%3D%202%20*%20half%20%2F%20npix%0A%20%20%20%20%20%20%20%20t0%20%3D%20time.perf_counter()%0A%20%20%20%20%20%20%20%20with%20rio.open(url)%20as%20probe%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20meta%20%3D%20dict(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20driver%3Dprobe.driver%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20size%3Df%22%7Bprobe.width%7D%20x%20%7Bprobe.height%7D%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20bands%3Dprobe.count%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20dtype%3Dprobe.dtypes%5B0%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20crs%3Dstr(probe.crs)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20overviews%3Dprobe.overviews(1)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20nodata%3Dprobe.nodata%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20ovr%20%3D%20_pick_overview(probe%2C%20res_m)%0A%20%20%20%20%20%20%20%20with%20rio.open(url%2C%20overview_level%3Dovr)%20as%20src%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20with%20WarpedVRT(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20src%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20crs%3Ddst_crs%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20transform%3Ddst_transform%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20width%3Dnpix%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20height%3Dnpix%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20resampling%3DResampling.bilinear%20if%20src.count%20%3E%201%20else%20Resampling.average%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20add_alpha%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%20as%20vrt%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20arr%20%3D%20vrt.read()%0A%20%20%20%20%20%20%20%20meta%5B%22overview_used%22%5D%20%3D%20ovr%0A%20%20%20%20%20%20%20%20meta%5B%22seconds%22%5D%20%3D%20round(time.perf_counter()%20-%20t0%2C%201)%0A%20%20%20%20%20%20%20%20meta%5B%22res_m%22%5D%20%3D%20res_m%0A%20%20%20%20%20%20%20%20return%20arr%2C%20meta%0A%0A%20%20%20%20warped%2C%20raster_meta%20%3D%20warp_square(SPILHAUS%2C%20HALF%2C%20NPIX)%0A%20%20%20%20dst_res_m%20%3D%20raster_meta%5B%22res_m%22%5D%0A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20f%22Warped%20%60%7BRASTER_URL.split('%2F')%5B-1%5D%7D%60%20to%20a%20%7BNPIX%7D%20x%20%7BNPIX%7D%20Spilhaus%20grid%20%22%0A%20%20%20%20%20%20%20%20f%22(%7Bdst_res_m%2F1000%3A.1f%7D%20km%2Fpx)%20in%20**%7Braster_meta%5B'seconds'%5D%7D%20s**%2C%20reading%20overview%20level%20%22%0A%20%20%20%20%20%20%20%20f%22%60%7Braster_meta%5B'overview_used'%5D%7D%60%20of%20%7Braster_meta%5B'overviews'%5D%20or%20'none'%7D.%22%0A%20%20%20%20)%0A%20%20%20%20return%20dst_res_m%2C%20rio%2C%20warp_square%2C%20warped%0A%0A%0A%40app.cell%0Adef%20_(HALF%2C%20dst_res_m%2C%20matplotlib%2C%20mo%2C%20np%2C%20plt%2C%20warped)%3A%0A%20%20%20%20from%20matplotlib.colors%20import%20LightSource%2C%20LinearSegmentedColormap%0A%0A%20%20%20%20def%20render_raster(ax%2C%20arr%2C%20half%2C%20res_m)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Draw%20the%20warped%20array.%20Single%20band%20-%3E%20bathymetry%20with%20land%20flat%3B%203%2B%20bands%20-%3E%20RGB.%22%22%22%0A%20%20%20%20%20%20%20%20alpha%20%3D%20arr%5B-1%5D%20%3E%200%0A%20%20%20%20%20%20%20%20if%20arr.shape%5B0%5D%20%3E%3D%204%3A%20%20%23%20RGB(%2Balpha)%0A%20%20%20%20%20%20%20%20%20%20%20%20rgb%20%3D%20np.moveaxis(arr%5B%3A3%5D%2C%200%2C%20-1).astype(float)%20%2F%20255.0%0A%20%20%20%20%20%20%20%20%20%20%20%20rgba%20%3D%20np.dstack(%5Brgb%2C%20alpha.astype(float)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20ax.imshow(rgba%2C%20extent%3D(-half%2C%20half%2C%20-half%2C%20half)%2C%20interpolation%3D%22nearest%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20%22rgb%22%0A%20%20%20%20%20%20%20%20z%20%3D%20arr%5B0%5D.astype(float)%0A%20%20%20%20%20%20%20%20z%5B~alpha%5D%20%3D%20np.nan%0A%20%20%20%20%20%20%20%20ocean%20%3D%20np.where(z%20%3C%200%2C%20z%2C%20np.nan)%0A%20%20%20%20%20%20%20%20land%20%3D%20np.where(z%20%3E%3D%200%2C%20z%2C%20np.nan)%0A%20%20%20%20%20%20%20%20%23%20ocean%3A%20deep-to-shallow%20blues%2C%20with%20a%20gentle%20hillshade%20for%20ridges%20and%20trenches%0A%20%20%20%20%20%20%20%20cmap%20%3D%20LinearSegmentedColormap.from_list(%0A%20%20%20%20%20%20%20%20%20%20%20%20%22bathy%22%2C%20%5B%22%2308143f%22%2C%20%22%23123c7a%22%2C%20%22%232a6fb5%22%2C%20%22%236aa9d8%22%2C%20%22%23b9d9ea%22%2C%20%22%23e3f1f6%22%5D%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20norm%20%3D%20matplotlib.colors.Normalize(vmin%3D-7000%2C%20vmax%3D0)%0A%20%20%20%20%20%20%20%20rgba%20%3D%20cmap(norm(np.nan_to_num(ocean%2C%20nan%3D0.0)))%0A%20%20%20%20%20%20%20%20ls%20%3D%20LightSource(azdeg%3D315%2C%20altdeg%3D45)%0A%20%20%20%20%20%20%20%20hs%20%3D%20ls.hillshade(np.nan_to_num(ocean%2C%20nan%3D0.0)%2C%20vert_exag%3D15.0%2C%20dx%3Dres_m%2C%20dy%3Dres_m)%0A%20%20%20%20%20%20%20%20rgba%5B...%2C%20%3A3%5D%20*%3D%20(0.65%20%2B%200.45%20*%20hs)%5B...%2C%20None%5D%0A%20%20%20%20%20%20%20%20rgba%5B...%2C%20%3A3%5D%20%3D%20rgba%5B...%2C%20%3A3%5D.clip(0%2C%201)%0A%20%20%20%20%20%20%20%20rgba%5B...%2C%203%5D%20%3D%20np.where(np.isfinite(ocean)%2C%201.0%2C%200.0)%0A%20%20%20%20%20%20%20%20ax.imshow(rgba%2C%20extent%3D(-half%2C%20half%2C%20-half%2C%20half)%2C%20interpolation%3D%22nearest%22)%0A%20%20%20%20%20%20%20%20land_rgba%20%3D%20np.zeros(z.shape%20%2B%20(4%2C))%0A%20%20%20%20%20%20%20%20land_rgba%5B...%2C%20%3A3%5D%20%3D%20matplotlib.colors.to_rgb(%22%23d9d2c5%22)%0A%20%20%20%20%20%20%20%20land_rgba%5B...%2C%203%5D%20%3D%20np.where(np.isfinite(land)%2C%201.0%2C%200.0)%0A%20%20%20%20%20%20%20%20ax.imshow(land_rgba%2C%20extent%3D(-half%2C%20half%2C%20-half%2C%20half)%2C%20interpolation%3D%22nearest%22)%0A%20%20%20%20%20%20%20%20return%20%22bathymetry%22%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%208))%0A%20%20%20%20raster_kind%20%3D%20render_raster(_ax%2C%20warped%2C%20HALF%2C%20dst_res_m)%0A%20%20%20%20_ax.set_xlim(-HALF%2C%20HALF)%0A%20%20%20%20_ax.set_ylim(-HALF%2C%20HALF)%0A%20%20%20%20_ax.set_axis_off()%0A%20%20%20%20_ax.set_title(f%22Raster%20only%20(%7Braster_kind%7D)%2C%20warped%20by%20GDAL%20via%20rasterio%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20mo.vstack(%5B_fig%5D)%0A%20%20%20%20return%20(render_raster%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%204.%20The%20same%20warp%20through%20other%20doors%0A%0A%20%20%20%20Same%20source%2C%20same%20target%20CRS%2C%20same%20grid%20-%20requested%20via%20each%20command-line%0A%20%20%20%20tool%20that%20happens%20to%20be%20on%20%60PATH%60.%20Success%20or%20failure%20here%20is%20*purely*%20a%0A%20%20%20%20function%20of%20which%20PROJ%20that%20binary%20links.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(HALF%2C%20RASTER_URL%2C%20SPILHAUS%2C%20mo%2C%20os%2C%20shutil%2C%20subprocess)%3A%0A%20%20%20%20import%20tempfile%0A%0A%20%20%20%20_out_dir%20%3D%20tempfile.mkdtemp(prefix%3D%22spilhaus_%22)%0A%20%20%20%20_small%20%3D%20256%20%20%23%20tiny%20output%3B%20this%20is%20a%20capability%20probe%2C%20not%20a%20render%0A%20%20%20%20_te%20%3D%20%5B%22-te%22%2C%20str(-HALF)%2C%20str(-HALF)%2C%20str(HALF)%2C%20str(HALF)%5D%0A%0A%20%20%20%20_cmds%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22gdalwarp%20(classic%20CLI)%22%3A%20%5B%22gdalwarp%22%2C%20%22-q%22%2C%20%22-overwrite%22%2C%20%22-t_srs%22%2C%20SPILHAUS%2C%20*_te%2C%20%22-ts%22%2C%20str(_small)%2C%20str(_small)%2C%20RASTER_URL%2C%20f%22%7B_out_dir%7D%2Fgdalwarp.tif%22%5D%2C%0A%20%20%20%20%20%20%20%20%22gdal%20raster%20reproject%20(unified%20CLI)%22%3A%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22gdal%22%2C%20%22raster%22%2C%20%22reproject%22%2C%20%22--overwrite%22%2C%20%22--dst-crs%22%2C%20SPILHAUS%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22--bbox%22%2C%20f%22%7B-HALF%7D%2C%7B-HALF%7D%2C%7BHALF%7D%2C%7BHALF%7D%22%2C%20%22--bbox-crs%22%2C%20SPILHAUS%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22--size%22%2C%20f%22%7B_small%7D%2C%7B_small%7D%22%2C%20RASTER_URL%2C%20f%22%7B_out_dir%7D%2Fgdal_cli.tif%22%2C%0A%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20%20%20%20%20%22rio%20warp%20(rasterio%20CLI)%22%3A%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22rio%22%2C%20%22warp%22%2C%20%22--overwrite%22%2C%20%22--dst-crs%22%2C%20SPILHAUS%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22--bounds%22%2C%20str(-HALF)%2C%20str(-HALF)%2C%20str(HALF)%2C%20str(HALF)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22--dimensions%22%2C%20str(_small)%2C%20str(_small)%2C%20RASTER_URL%2C%20f%22%7B_out_dir%7D%2Frio.tif%22%2C%0A%20%20%20%20%20%20%20%20%5D%2C%0A%20%20%20%20%7D%0A%0A%20%20%20%20cli_results%20%3D%20%5B%5D%0A%20%20%20%20for%20_label%2C%20_cmd%20in%20_cmds.items()%3A%0A%20%20%20%20%20%20%20%20_exe%20%3D%20shutil.which(_cmd%5B0%5D)%0A%20%20%20%20%20%20%20%20if%20not%20_exe%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20cli_results.append(dict(tool%3D_label%2C%20result%3D%22not%20on%20PATH%22%2C%20detail%3D%22%22))%0A%20%20%20%20%20%20%20%20%20%20%20%20continue%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_p%20%3D%20subprocess.run(_cmd%2C%20capture_output%3DTrue%2C%20text%3DTrue%2C%20timeout%3D600)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20_p.returncode%20%3D%3D%200%20and%20os.path.exists(_cmd%5B-1%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cli_results.append(dict(tool%3D_label%2C%20result%3D%22ok%22%2C%20detail%3Df%22%7Bos.path.getsize(_cmd%5B-1%5D)%2F%2F1024%7D%20KiB%20written%22))%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_err%20%3D%20(_p.stderr%20or%20_p.stdout).strip().splitlines()%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20cli_results.append(dict(tool%3D_label%2C%20result%3D%22FAIL%22%2C%20detail%3D_err%5B-1%5D%20if%20_err%20else%20f%22exit%20%7B_p.returncode%7D%22))%0A%20%20%20%20%20%20%20%20except%20subprocess.TimeoutExpired%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20cli_results.append(dict(tool%3D_label%2C%20result%3D%22timeout%22%2C%20detail%3D%22%22))%0A%0A%20%20%20%20mo.ui.table(cli_results%2C%20selection%3DNone%2C%20label%3D%22Command-line%20warp%20probes%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%205.%20Vector%3A%20Natural%20Earth%20ocean%20polygons%0A%0A%20%20%20%20Read%20remotely%20with%20%60pyogrio%60%20(via%20%60geopandas.read_file%60)%20straight%20from%20the%0A%20%20%20%20Natural%20Earth%20GitHub%20repository%20with%20%60%2Fvsicurl%2F%60%20-%20no%20download%20step.%0A%0A%20%20%20%20Why%20*ocean*%20rather%20than%20*land*%3A%20the%20Spilhaus%20cut%20runs%20through%20the%0A%20%20%20%20continents%2C%20so%20ocean%20rings%20barely%20straddle%20it.%20The%20few%20segments%20that%20do%0A%20%20%20%20cross%20are%20detected%20as%20absurdly%20long%20jumps%20in%20projected%20space%20and%20dropped%2C%0A%20%20%20%20and%20coastlines%20are%20drawn%20as%20lines%20rather%20than%20filled%20polygons.%20(A%20polygon%0A%20%20%20%20fill%20would%20need%20the%20rings%20split%20along%20the%20cut%20first%3B%20that%20is%20a%20separate%2C%0A%20%20%20%20interesting%20exercise.)%0A%0A%20%20%20%20%60geopandas.to_crs%60%20delegates%20to%20%60pyproj%60%3B%20if%20that%20PROJ%20predates%209.6%2C%20it%0A%20%20%20%20cannot%20help%20here%2C%20so%20the%20transform%20goes%20coordinate-by-coordinate%20through%0A%20%20%20%20whichever%20engine%20section%201%20found%2C%20using%20%60shapely.transform%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(OCEAN_URL%2C%20forward%2C%20np)%3A%0A%20%20%20%20import%20geopandas%20as%20gpd%0A%20%20%20%20import%20shapely%0A%0A%20%20%20%20ocean%20%3D%20gpd.read_file(OCEAN_URL)%0A%0A%20%20%20%20def%20rings_to_segments(gdf%2C%20densify_deg%3D0.5%2C%20max_jump%3DNone%2C%20fwd%3Dforward)%3A%0A%20%20%20%20%20%20%20%20%22%22%22All%20polygon%20rings%20-%3E%20list%20of%20(N%2C2)%20projected%20coordinate%20arrays%2C%20split%20at%20long%20jumps.%22%22%22%0A%20%20%20%20%20%20%20%20out%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20geom%20in%20gdf.geometry%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20geom%20%3D%20shapely.segmentize(geom%2C%20densify_deg)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20poly%20in%20getattr(geom%2C%20%22geoms%22%2C%20%5Bgeom%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20ring%20in%20%5Bpoly.exterior%2C%20*poly.interiors%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20c%20%3D%20np.asarray(ring.coords)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20x%2C%20y%20%3D%20fwd(c%5B%3A%2C%200%5D%2C%20c%5B%3A%2C%201%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20xy%20%3D%20np.c_%5Bx%2C%20y%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20step%20%3D%20np.hypot(np.diff(x)%2C%20np.diff(y))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20segments%20that%20are%20artefacts%20of%20the%20lon%2Flat%20representation%2C%20not%20coastline%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20both%20ends%20on%20the%20antimeridian%2C%20or%20both%20ends%20at%20a%20pole%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20lon0%2C%20lon1%20%3D%20c%5B%3A-1%2C%200%5D%2C%20c%5B1%3A%2C%200%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20lat0%2C%20lat1%20%3D%20c%5B%3A-1%2C%201%5D%2C%20c%5B1%3A%2C%201%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20synthetic%20%3D%20((np.abs(lon0)%20%3E%20179.99)%20%26%20(np.abs(lon1)%20%3E%20179.99))%20%7C%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(np.abs(lat0)%20%3E%2089.99)%20%26%20(np.abs(lat1)%20%3E%2089.99)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ok%20%3D%20np.isfinite(step)%20%26%20(step%20%3C%20max_jump)%20%26%20~synthetic%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20split%20into%20runs%20of%20consecutive%20good%20segments%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20idx%20%3D%20np.flatnonzero(~ok)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20start%20%3D%200%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20i%20in%20list(idx)%20%2B%20%5Blen(step)%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20i%20-%20start%20%3E%3D%201%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20out.append(xy%5Bstart%20%3A%20i%20%2B%201%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20start%20%3D%20i%20%2B%201%0A%20%20%20%20%20%20%20%20return%20out%0A%0A%20%20%20%20return%20ocean%2C%20rings_to_segments%0A%0A%0A%40app.cell%0Adef%20_(ENGINE%2C%20HALF%2C%20mo%2C%20ocean%2C%20rings_to_segments)%3A%0A%20%20%20%20coast_segments%20%3D%20rings_to_segments(ocean%2C%20densify_deg%3D0.5%2C%20max_jump%3DHALF%20%2F%208)%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20f%22Ocean%20layer%3A%20%7Blen(ocean)%7D%20feature(s)%2C%20CRS%20%60%7Bocean.crs%7D%60%3B%20%22%0A%20%20%20%20%20%20%20%20f%22%7Blen(coast_segments)%7D%20coastline%20runs%20after%20projecting%20with%20**%7BENGINE%7D**%20and%20dropping%20cut-crossing%20segments.%22%0A%20%20%20%20)%0A%20%20%20%20return%20(coast_segments%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%206.%20Points%3A%20landmarks%2C%20transformed%20directly%0A%0A%20%20%20%20No%20GeoDataFrame%2C%20no%20raster%3A%20a%20dict%20of%20lon%2Flat%20pairs%20pushed%20through%20the%0A%20%20%20%20forward%20transform.%20This%20is%20the%20smallest%20possible%20%22does%20my%20PROJ%20have%0A%20%20%20%20Spilhaus%22%20test%2C%20and%20the%20one%20to%20paste%20into%20a%20bug%20report.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(forward%2C%20mo)%3A%0A%20%20%20%20landmarks%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22Hobart%22%3A%20(147.33%2C%20-42.88)%2C%0A%20%20%20%20%20%20%20%20%22Tokyo%22%3A%20(139.69%2C%2035.69)%2C%0A%20%20%20%20%20%20%20%20%22Reykjavik%22%3A%20(-21.94%2C%2064.15)%2C%0A%20%20%20%20%20%20%20%20%22Cape%20of%20Good%20Hope%22%3A%20(18.47%2C%20-34.36)%2C%0A%20%20%20%20%20%20%20%20%22Bering%20Strait%22%3A%20(-169.0%2C%2065.8)%2C%0A%20%20%20%20%20%20%20%20%22Drake%20Passage%22%3A%20(-63.0%2C%20-59.0)%2C%0A%20%20%20%20%20%20%20%20%22Point%20Nemo%22%3A%20(-123.39%2C%20-48.88)%2C%0A%20%20%20%20%20%20%20%20%22Challenger%20Deep%22%3A%20(142.20%2C%2011.37)%2C%0A%20%20%20%20%20%20%20%20%22Puerto%20Rico%20Trench%22%3A%20(-66.5%2C%2019.7)%2C%0A%20%20%20%20%20%20%20%20%22Mid-Atlantic%20Ridge%20(Azores)%22%3A%20(-27.0%2C%2038.0)%2C%0A%20%20%20%20%20%20%20%20%22Macquarie%20Island%22%3A%20(158.94%2C%20-54.62)%2C%0A%20%20%20%20%20%20%20%20%22Mawson%20Station%22%3A%20(62.87%2C%20-67.60)%2C%0A%20%20%20%20%20%20%20%20%22Casey%20Station%22%3A%20(110.53%2C%20-66.28)%2C%0A%20%20%20%20%7D%0A%20%20%20%20_lon%20%3D%20%5Bv%5B0%5D%20for%20v%20in%20landmarks.values()%5D%0A%20%20%20%20_lat%20%3D%20%5Bv%5B1%5D%20for%20v%20in%20landmarks.values()%5D%0A%20%20%20%20_x%2C%20_y%20%3D%20forward(_lon%2C%20_lat)%0A%20%20%20%20landmark_xy%20%3D%20%7Bk%3A%20(float(x)%2C%20float(y))%20for%20k%2C%20x%2C%20y%20in%20zip(landmarks%2C%20_x%2C%20_y)%7D%0A%20%20%20%20_rows%20%3D%20%5Bdict(name%3Dk%2C%20lon%3Dv%5B0%5D%2C%20lat%3Dv%5B1%5D%2C%20x%3Dround(landmark_xy%5Bk%5D%5B0%5D)%2C%20y%3Dround(landmark_xy%5Bk%5D%5B1%5D))%20for%20k%2C%20v%20in%20landmarks.items()%5D%0A%20%20%20%20mo.ui.table(_rows%2C%20selection%3DNone%2C%20label%3D%22Landmarks%20in%20Spilhaus%20metres%22)%0A%20%20%20%20return%20landmark_xy%2C%20landmarks%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%207.%20Composite%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20ENGINE%2C%0A%20%20%20%20HALF%2C%0A%20%20%20%20coast_segments%2C%0A%20%20%20%20dst_res_m%2C%0A%20%20%20%20landmark_xy%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20probes%2C%0A%20%20%20%20render_raster%2C%0A%20%20%20%20warped%2C%0A)%3A%0A%20%20%20%20from%20matplotlib%20import%20patheffects%0A%20%20%20%20from%20matplotlib.collections%20import%20LineCollection%0A%0A%20%20%20%20fig%2C%20ax%20%3D%20plt.subplots(figsize%3D(11%2C%2011))%0A%20%20%20%20ax.set_facecolor(%22%23d9d2c5%22)%0A%20%20%20%20render_raster(ax%2C%20warped%2C%20HALF%2C%20dst_res_m)%0A%20%20%20%20ax.add_collection(LineCollection(coast_segments%2C%20colors%3D%22%233a3a3a%22%2C%20linewidths%3D0.45))%0A%20%20%20%20for%20_name%2C%20(_x%2C%20_y)%20in%20landmark_xy.items()%3A%0A%20%20%20%20%20%20%20%20ax.plot(_x%2C%20_y%2C%20%22o%22%2C%20ms%3D5%2C%20color%3D%22%23f2b134%22%2C%20mec%3D%22black%22%2C%20mew%3D0.6%2C%20zorder%3D5)%0A%20%20%20%20%20%20%20%20ax.annotate(%0A%20%20%20%20%20%20%20%20%20%20%20%20_name%2C%20(_x%2C%20_y)%2C%20xytext%3D(5%2C%204)%2C%20textcoords%3D%22offset%20points%22%2C%20fontsize%3D8.5%2C%20zorder%3D6%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20annotation_clip%3DTrue%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20path_effects%3D%5Bpatheffects.withStroke(linewidth%3D2.5%2C%20foreground%3D%22white%22)%5D%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20ax.set_xlim(-HALF%2C%20HALF)%0A%20%20%20%20ax.set_ylim(-HALF%2C%20HALF)%0A%20%20%20%20ax.set_axis_off()%0A%20%20%20%20_proj_used%20%3D%20next((p%5B%22proj%22%5D%20for%20p%20in%20probes%20if%20p%5B%22package%22%5D%20%3D%3D%20ENGINE)%2C%20%22%3F%22)%0A%20%20%20%20ax.set_title(%0A%20%20%20%20%20%20%20%20f%22The%20world%20ocean%20in%20Spilhaus%20projection%5Cn%22%0A%20%20%20%20%20%20%20%20f%22raster%20warped%20by%20GDAL%20(rasterio)%3B%20coastlines%20and%20points%20transformed%20by%20%7BENGINE%7D%20(PROJ%20%7B_proj_used%7D)%22%2C%0A%20%20%20%20%20%20%20%20fontsize%3D12%2C%0A%20%20%20%20%20%20%20%20pad%3D16%2C%0A%20%20%20%20)%0A%20%20%20%20fig.tight_layout()%0A%20%20%20%20fig.savefig(%22spilhaus_composite.png%22%2C%20dpi%3D150%2C%20facecolor%3D%22white%22)%0A%20%20%20%20mo.vstack(%5Bfig%2C%20mo.md(%22Saved%20as%20%60spilhaus_composite.png%60%20next%20to%20the%20notebook.%22)%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%208.%20Naming%20the%20thing%3A%20proj-string%2C%20WKT%2C%20and%20the%20missing%20authority%20code%0A%0A%20%20%20%20%60%2Bproj%3Dspilhaus%60%20is%20a%20fine%20thing%20to%20type.%20It%20is%20a%20poor%20thing%20to%20store%2C%20cite%2C%0A%20%20%20%20or%20hand%20to%20someone%20in%20five%20years.%20Three%20renderings%20of%20%22the%20same%22%20projection%3A%0A%0A%20%20%20%201.%20the%20bare%20proj-string%2C%20exported%20to%20WKT2%20-%20the%20method%20is%20named%2C%20but%20every%0A%20%20%20%20%20%20%20parameter%20is%20*implicit*%3A%20the%20numbers%20live%20in%20PROJ's%20source%20code%3B%0A%20%20%20%202.%20the%20same%20proj-string%20with%20all%20parameters%20spelled%20out%20-%20now%20the%20WKT%20carries%0A%20%20%20%20%20%20%20the%20numbers%2C%20and%20round-trips%3B%0A%20%20%20%203.%20**ESRI%3A54099**%2C%20the%20only%20authority%20code%20you%20will%20find.%20It%20is%20a%20*different*%20map%3A%0A%20%20%20%20%20%20%20ESRI's%20Adams%20Square%20II%20convention%20has%20the%20square%20a%20factor%20of%20sqrt(2)%20larger%2C%0A%20%20%20%20%20%20%20which%20PROJ%20expresses%20by%20mapping%20it%20to%20%60%2Bk_0%3D1.41421356237%60.%0A%0A%20%20%20%20There%20is%20no%20EPSG%20code.%20Nobody%20has%20registered%20one%2C%20and%20for%20a%20projection%20whose%0A%20%20%20%20whole%20point%20is%20one%20specific%20choice%20of%20parameters%20that%20is%20a%20gap%20worth%20noticing.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(SPILHAUS%2C%20SPILHAUS_ALT%2C%20SPILHAUS_EXPLICIT%2C%20mo%2C%20rio)%3A%0A%20%20%20%20from%20rasterio.warp%20import%20transform%20as%20_xf%0A%0A%20%20%20%20def%20_wkt(crs_in)%3A%0A%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20crs%20%3D%20rio.crs.CRS.from_string(crs_in)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20crs%2C%20crs.to_wkt(version%3D%22WKT2_2019%22%2C%20pretty%3DTrue)%0A%20%20%20%20%20%20%20%20except%20Exception%20as%20e%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20None%2C%20f%22(%7Btype(e).__name__%7D%3A%20%7Be%7D)%22%0A%0A%20%20%20%20_hobart%20%3D%20(147.33%2C%20-42.88)%0A%20%20%20%20renderings%20%3D%20%5B%5D%0A%20%20%20%20for%20_label%2C%20_src%20in%20%5B%0A%20%20%20%20%20%20%20%20(%22bare%20proj-string%22%2C%20SPILHAUS)%2C%0A%20%20%20%20%20%20%20%20(%22explicit%20proj-string%22%2C%20SPILHAUS_EXPLICIT)%2C%0A%20%20%20%20%20%20%20%20(%22ESRI%3A54099%22%2C%20%22ESRI%3A54099%22)%2C%0A%20%20%20%20%20%20%20%20(%22alternative%20centre%20(section%209)%22%2C%20SPILHAUS_ALT)%2C%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20_crs%2C%20_w%20%3D%20_wkt(_src)%0A%20%20%20%20%20%20%20%20_xy%20%3D%20%22-%22%0A%20%20%20%20%20%20%20%20_p4%20%3D%20%22-%22%0A%20%20%20%20%20%20%20%20if%20_crs%20is%20not%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_x%2C%20_y%20%3D%20_xf(%22EPSG%3A4326%22%2C%20_crs%2C%20%5B_hobart%5B0%5D%5D%2C%20%5B_hobart%5B1%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_xy%20%3D%20f%22%7B_x%5B0%5D%3A%2C.0f%7D%2C%20%7B_y%5B0%5D%3A%2C.0f%7D%22%0A%20%20%20%20%20%20%20%20%20%20%20%20except%20Exception%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_xy%20%3D%20%22transform%20failed%22%0A%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_p4%20%3D%20_crs.to_proj4()%0A%20%20%20%20%20%20%20%20%20%20%20%20except%20Exception%3A%20%20%23%20noqa%3A%20BLE001%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_p4%20%3D%20%22(no%20proj-string%20export)%22%0A%20%20%20%20%20%20%20%20renderings.append(dict(rendering%3D_label%2C%20input%3D_src%2C%20proj_string_export%3D_p4%2C%20hobart_xy%3D_xy%2C%20wkt2%3D_w))%0A%0A%20%20%20%20mo.vstack(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.ui.table(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%5B%7Bk%3A%20v%20for%20k%2C%20v%20in%20r.items()%20if%20k%20!%3D%20%22wkt2%22%7D%20for%20r%20in%20renderings%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20selection%3DNone%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22Same%20projection%2C%20several%20names.%20hobart_xy%20is%20Hobart's%20projected%20coordinate%20under%20each.%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.accordion(%7Bf%22WKT2%3A2019%20for%20%7Br%5B'rendering'%5D%7D%22%3A%20mo.md(f%22%60%60%60%5Cn%7Br%5B'wkt2'%5D%7D%5Cn%60%60%60%22)%20for%20r%20in%20renderings%7D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Look%20for%20%60PARAMETER%5B...%5D%60%20entries%20in%20the%20first%20WKT%3A%20there%20are%20none.%20The%20map%20is%20fully%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22determined%20only%20by%20defaults%20in%20PROJ.%20The%20second%20rendering%20is%20what%20a%20technical%20citation%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22should%20carry.%20ESRI%3A54099%20places%20Hobart%20somewhere%20else%20entirely%3B%20it%20is%20a%20bigger%20square.%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20kind%3D%22info%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%209.%20Another%20member%20of%20the%20family%0A%0A%20%20%20%20Spilhaus's%20genius%20was%20not%20the%20projection%20method%20(it%20is%20Adams'%20World%20in%20a%20Square%20II%2C%0A%20%20%20%201929)%20but%20the%20*parameters*%3A%20a%20centre%2C%20an%20azimuth%20and%20a%20rotation%20chosen%20so%20the%20cut%20falls%0A%20%20%20%20almost%20entirely%20on%20land.%20Change%20them%20and%20you%20have%20a%20perfectly%20valid%20map%20of%20the%20same%0A%20%20%20%20family%20that%20nobody%20would%20call%20the%20Spilhaus%20map.%20Below%2C%20a%20naive%20recentring%20on%20the%0A%20%20%20%20Atlantic.%20Same%20method%2C%20same%20code%20path%2C%20and%20the%20square's%20edges%20now%20slice%20the%20ocean.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20HALF%2C%0A%20%20%20%20NPIX%2C%0A%20%20%20%20SPILHAUS_ALT%2C%0A%20%20%20%20landmarks%2C%0A%20%20%20%20make_forward%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20ocean%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20render_raster%2C%0A%20%20%20%20rings_to_segments%2C%0A%20%20%20%20square_half%2C%0A%20%20%20%20warp_square%2C%0A)%3A%0A%20%20%20%20from%20matplotlib%20import%20patheffects%20as%20_pe%0A%20%20%20%20from%20matplotlib.collections%20import%20LineCollection%20as%20_LC%0A%0A%20%20%20%20forward_alt%20%3D%20make_forward(SPILHAUS_ALT)%0A%20%20%20%20HALF_ALT%20%3D%20square_half(forward_alt)%0A%20%20%20%20_npix%20%3D%20max(400%2C%20NPIX%20%2F%2F%202)%0A%20%20%20%20warped_alt%2C%20meta_alt%20%3D%20warp_square(SPILHAUS_ALT%2C%20HALF_ALT%2C%20_npix)%0A%20%20%20%20coast_alt%20%3D%20rings_to_segments(ocean%2C%20densify_deg%3D0.5%2C%20max_jump%3DHALF_ALT%20%2F%208%2C%20fwd%3Dforward_alt)%0A%20%20%20%20_lx%2C%20_ly%20%3D%20forward_alt(%5Bv%5B0%5D%20for%20v%20in%20landmarks.values()%5D%2C%20%5Bv%5B1%5D%20for%20v%20in%20landmarks.values()%5D)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%208))%0A%20%20%20%20_ax.set_facecolor(%22%23d9d2c5%22)%0A%20%20%20%20render_raster(_ax%2C%20warped_alt%2C%20HALF_ALT%2C%20meta_alt%5B%22res_m%22%5D)%0A%20%20%20%20_ax.add_collection(_LC(coast_alt%2C%20colors%3D%22%233a3a3a%22%2C%20linewidths%3D0.45))%0A%20%20%20%20_ax.plot(_lx%2C%20_ly%2C%20%22o%22%2C%20ms%3D4%2C%20color%3D%22%23f2b134%22%2C%20mec%3D%22black%22%2C%20mew%3D0.6%2C%20zorder%3D5)%0A%20%20%20%20for%20_name%2C%20_x%2C%20_y%20in%20zip(landmarks%2C%20_lx%2C%20_ly)%3A%0A%20%20%20%20%20%20%20%20_ax.annotate(%0A%20%20%20%20%20%20%20%20%20%20%20%20_name%2C%20(_x%2C%20_y)%2C%20xytext%3D(4%2C%203)%2C%20textcoords%3D%22offset%20points%22%2C%20fontsize%3D7.5%2C%20zorder%3D6%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20annotation_clip%3DTrue%2C%20path_effects%3D%5B_pe.withStroke(linewidth%3D2%2C%20foreground%3D%22white%22)%5D%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20_ax.set_xlim(-HALF_ALT%2C%20HALF_ALT)%0A%20%20%20%20_ax.set_ylim(-HALF_ALT%2C%20HALF_ALT)%0A%20%20%20%20_ax.set_axis_off()%0A%20%20%20%20_ax.set_title(f%22%7BSPILHAUS_ALT%7D%5Cn(half-width%20%7BHALF_ALT%2F1e6%3A.3f%7D%20Mm%20vs%20%7BHALF%2F1e6%3A.3f%7D%20Mm%20for%20the%20default)%22%2C%20fontsize%3D10%2C%20pad%3D12)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20mo.vstack(%5B_fig%5D)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%2010.%20What%20to%20take%20from%20this%0A%0A%20%20%20%20-%20A%20projection%20%22existing%20in%20PROJ%22%20is%20necessary%2C%20not%20sufficient.%20Each%20package%0A%20%20%20%20%20%20that%20bundles%20PROJ%20picks%20up%20the%20feature%20only%20when%20*its*%20wheel%20is%20rebuilt%20against%0A%20%20%20%20%20%20a%20new%20enough%20PROJ%2C%20on%20its%20own%20release%20cadence.%20The%20provenance%20table%20in%20section%201%20is%0A%20%20%20%20%20%20the%20actual%20answer%20to%20%22can%20I%20use%20Spilhaus%20here%22.%0A%20%20%20%20-%20GDAL-based%20paths%20(%60rasterio%60%2C%20%60pyogrio%60%2C%20the%20%60osgeo%60%20bindings%2C%20CLI%20tools)%20each%0A%20%20%20%20%20%20link%20a%20PROJ%20too%3B%20they%20do%20not%20share%20one%20with%20%60pyproj%60.%0A%20%20%20%20-%20On%20conda-forge%20all%20of%20these%20share%20a%20single%20%60libproj%60%2C%20so%20the%20table%20collapses%20to%20one%0A%20%20%20%20%20%20row%20of%20truth.%20That%20is%20a%20real%20difference%20between%20the%20two%20ecosystems%2C%20and%20worth%20knowing%0A%20%20%20%20%20%20when%20someone%20asks%20why%20a%20projection%20%22works%20in%20QGIS%20but%20not%20in%20my%20notebook%22.%0A%20%20%20%20-%20The%20pattern%20generalises%3A%20any%20new%20PROJ%20operation%2C%20any%20new%20GDAL%20driver%20or%20CLI%20subcommand%0A%20%20%20%20%20%20(%60gdal%20raster%20reproject%60%2C%20section%204)%20flows%20downstream%20on%20the%20same%20uneven%20schedule.%0A%20%20%20%20-%20A%20proj-string%20is%20a%20recipe%20with%20hidden%20defaults.%20When%20it%20matters%2C%20write%20the%20WKT%20with%20the%0A%20%20%20%20%20%20parameters%20in%20it%20(section%208)%2C%20and%20do%20not%20assume%20the%20one%20authority%20code%20you%20can%20find%0A%20%20%20%20%20%20describes%20the%20map%20you%20have%20in%20mind.%0A%0A%20%20%20%20Re-run%20this%20notebook%20after%20every%20%60pip%20install%20-U%60%20and%20watch%20the%20rows%20change.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A