{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Applying the Ocean function to geometries and generating grids\n", "\n", "A common application is to test whether a certain geometry lies over the ocean and/or to generate an rasterized ocean mask.\n", "\n", "The [natural earth collection](https://www.naturalearthdata.com) contains ocean geometries. However, naively testing for inclusion of a geometry can take a very long time since the polygon associated with Asia a has many points and is not optimized for inclusion tests. Furthermore, for most applications, the isolated Capsian Sea may also not be considered to belong to the ocean.\n", "\n", "For these reasons a dedicated table can be generated in geoslurp, which will speed up inclusion searches and raster generation. This will be described below.\n", "\n", "The idea is that the coastal polygons are subdivided in smaller manageable pieces with PostGIS ([ST_subdivide](https://postgis.net/docs/ST_Subdivide.html)) and a geospatial index is generated to speed up queries." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", "%autoreload 2" ] }, { "cell_type": "code", "execution_count": 131, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Geoslurp-INFO: Reading CF convention defaults from /home/roelof/.cf-conventions.yaml\n" ] } ], "source": [ "from geoslurp.config.catalogue import DatasetCatalogue\n", "from geoslurp.db import Settings\n", "from geoslurp.config import setInfoLevel\n", "\n", "from geoslurp.db import geoslurpConnect\n", "from geoslurp.tools.shapelytools import gdal2rastio\n", "from geoslurp.tools.pandas import *\n", "from geoslurp.tools.cf import *\n", "from math import floor\n", "import xarray as xr\n", "import numpy as np\n", "setInfoLevel()\n", "\n", "gpcon=geoslurpConnect(readonly_user=False) # this will be a connection based on the readonly userfrom geoslurp.db.geo\n", "\n", "#Some datasets need info from the server side settings so we need to load these\n", "conf=Settings(gpcon)\n", "\n", "#catalogue\n", "gscat=DatasetCatalogue()\n", "\n", "tocea=\"natearth.ne_10m_oceanfunc\"\n", "#original table\n", "tocea_orig=\"natearth.ne_10m_ocean\"\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's pull and register the dedicated ocean function table\n", "```\n", "geoslurper.py -v --pull --register \"natearth.ne_10m_oceanfunc\"\n", "```\n", "\n", "Or alternatively in Python:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Geoslurp-INFO: using cached github catalogue /home/data/geoslurp_cache/githubcache/strawpants_geoshapes.yaml\n", "Geoslurp-INFO: using cached github catalogue /home/data/geoslurp_cache/githubcache/nvkelso_natural-earth-vector.yaml\n", "/usr/lib/python3.12/site-packages/osgeo/osr.py:410: FutureWarning: Neither osr.UseExceptions() nor osr.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default.\n", " warnings.warn(\n", "Geoslurp-INFO: creating table natearth.ne_10m_oceanfunc and index, this can take a while..\n", "Geoslurp-INFO: Done..\n" ] } ], "source": [ "dsoceanf=gscat.getDsetClass(conf, tocea)(gpcon)\n", "dsoceanf.pull()\n", "dsoceanf.register()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1: Generate a global 0.125 degree equidistant raster\n", "\n", "Note: The generated grid can also be found on [this github repository](https://github.com/strawpants/geoshapes)" ] }, { "cell_type": "code", "execution_count": 128, "metadata": {}, "outputs": [], "source": [ "resolution=0.125\n", "\n", "east=180\n", "west=-180\n", "south=-90\n", "north=90\n", "width=floor((east-west)/resolution)\n", "height=floor((north-south)/resolution)\n", "\n", "# create a query which creates a reference raster object and only selects those pixel longitude/lattiudes which lie within the polygons\n", "qry=f\"\"\"\n", " CREATE TEMPORARY TABLE oceanrast AS \n", " (SELECT ST_addband(ST_MakeEmptyRaster({width},{height},{west},{north},{resolution},-{resolution},0,0,srid=>4326),'8BUI'::text,1.0) as rast);\n", " SELECT ST_x(pnts.geom) as longitude ,ST_y(pnts.geom) as latitude from {tocea} as ocean INNER JOIN (SELECT (ST_PixelAsCentroids(rast, 1)).* from oceanrast ) pnts ON ST_Within(pnts.geom,ocean.geom);\n", "\"\"\"\n", "\n", "#load the valid lon latitudes in a pandas dataframe\n", "dfocean=pd.DataFrame.gslrp.load(gpcon,qry)\n", "\n", "#create an zero values grid in xaray\n", "lat=np.arange(south+resolution/2,north,resolution)\n", "lon=np.arange(west+resolution/2,east,resolution)\n", "\n", "daocean=xr.DataArray(data=0.0,dims=[\"latitude\",\"longitude\"],coords={\"longitude\":(\"longitude\",lon),\"latitude\":(\"latitude\",lat)}).stack(lonlat=[\"longitude\",\"latitude\"])\n", "#assign 1 at ocean points\n", "daocean.loc[pd.MultiIndex.from_frame(dfocean)]=1.0\n", "daocean=daocean.unstack('lonlat').T" ] }, { "cell_type": "code", "execution_count": 143, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.Dataset>\n",
"Dimensions: (longitude: 2880, latitude: 1440)\n",
"Coordinates:\n",
" * longitude (longitude) float64 -179.9 -179.8 -179.7 ... 179.7 179.8 179.9\n",
" * latitude (latitude) float64 -89.94 -89.81 -89.69 ... 89.69 89.81 89.94\n",
" spatial_ref int64 0\n",
"Data variables:\n",
" oceanfunc (latitude, longitude) float64 0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0\n",
"Attributes:\n",
" Conventions: CF-1.9\n",
" title: Gridded Ocean mask generated from natural earth 10m dataset\n",
" institution: Roelof Rietbroek <r.rietbroek@utwente.nl>, Faculty of Geoin...\n",
" source: geoslurp\n",
" history: 2024-11-26 14:53:20.439413 geoslurp\n",
" references: https://geoslurp.wobbly.earth/en/latest/notebooks/OceanFunc...\n",
" comment: Auto generated