{ "cells": [ { "cell_type": "markdown", "id": "9fa1910a-6db7-4592-986c-e80ac015ac26", "metadata": {}, "source": [ "# Calculate distances and find neighbors" ] }, { "cell_type": "code", "execution_count": 1, "id": "a293cd46-83e6-4d73-adcc-ee9a6724739f", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "\n", "import poligrain as plg" ] }, { "cell_type": "markdown", "id": "ee6c19b3-d6a9-4669-b7c8-53894c1ceee6", "metadata": {}, "source": [ "## Point-to-point distances" ] }, { "cell_type": "markdown", "id": "19f3e1b8-51ad-4a72-8413-7cb619c9dee3", "metadata": {}, "source": [ "### Get PWS dataset" ] }, { "cell_type": "code", "execution_count": 2, "id": "df67e076-5ca4-43b2-92dc-552d19efc32a", "metadata": {}, "outputs": [], "source": [ "!curl -OL https://github.com/OpenSenseAction/training_school_opensene_2023/raw/main/data/pws/data_PWS_netCDF_AMS_float.nc" ] }, { "cell_type": "code", "execution_count": 3, "id": "61b86248-8c91-4a9d-ba26-a3e168d31945", "metadata": {}, "outputs": [], "source": [ "ds_pws = xr.open_dataset(\"data_PWS_netCDF_AMS_float.nc\")\n", "\n", "# fix some issues with this dataset\n", "ds_pws[\"time\"] = pd.to_datetime(ds_pws.time.data, unit=\"s\")\n", "ds_pws[\"lon\"] = (\"id\", ds_pws.lon.data)\n", "ds_pws[\"lat\"] = (\"id\", ds_pws.lat.data)" ] }, { "cell_type": "markdown", "id": "bf4ef4f3-640a-4c9c-b32f-8af551ee0821", "metadata": {}, "source": [ "### Project cooridnates from lon-lat to UTM zone for Europe\n", "\n", "To do meaningful distance calculations we have to project the lon-lat coordinates first.\n", "\n", "Info on the UTM zone 32N projection `'EPSG:25832'` that is used can be found [here](https://epsg.io/25832#google_vignette)\n", "\n", "Note that `plg.spatial.project_coordinates` can also use a different source coordinate system, but its default is WGS 84 (lon and lat in degrees), which is `'EPSG:4326'` " ] }, { "cell_type": "code", "execution_count": 4, "id": "5040ff46-0608-4496-8e08-60ceb42c3464", "metadata": {}, "outputs": [], "source": [ "ds_pws.coords[\"x\"], ds_pws.coords[\"y\"] = plg.spatial.project_point_coordinates(\n", " ds_pws.lon, ds_pws.lat, \"EPSG:25832\"\n", ")" ] }, { "cell_type": "code", "execution_count": 5, "id": "89ea4f3f-e2c2-46ff-885c-e518170ef69b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.Dataset> Size: 237MB\n",
"Dimensions: (time: 219168, id: 134)\n",
"Coordinates:\n",
" * time (time) datetime64[ns] 2MB 2016-05-01T00:05:00 ... 2018-06-01\n",
" * id (id) <U6 3kB 'ams1' 'ams2' 'ams3' ... 'ams132' 'ams133' 'ams134'\n",
" elevation (id) <U3 2kB ...\n",
" lat (id) float64 1kB 52.31 52.3 52.31 52.35 ... 52.43 52.3 52.26\n",
" lon (id) float64 1kB 4.671 4.675 4.677 4.678 ... 5.036 5.041 5.045\n",
" x (id) float64 1kB 2.049e+05 2.052e+05 ... 2.301e+05 2.301e+05\n",
" y (id) float64 1kB 5.804e+06 5.803e+06 ... 5.802e+06 5.798e+06\n",
"Data variables:\n",
" rainfall (id, time) float64 235MB ...\n",
"Attributes:\n",
" title: PWS data from Amsterdam\n",
" institution: Wageningen University and Research, Department of Environm...\n",
" history: Test Version 0.1\n",
" references: https://doi.org/10.1029/2019GL083731\n",
" date_created: 2022-10-18 10:32:00\n",
" Conventions: OPENSENSE V0\n",
" location: Amsterdam Metropolitan Area\n",
" source: Netamo\n",
" comment: <xarray.DataArray (id: 134, id_neighbor: 134)> Size: 144kB\n",
"array([[ 0. , 518.79828487, 531.99603941, ...,\n",
" 28532.23134009, 25291.03668667, 25994.5064412 ],\n",
" [ 518.79828487, 0. , 728.38027968, ...,\n",
" 28493.8273213 , 24995.94592241, 25633.14869253],\n",
" [ 531.99603941, 728.38027968, 0. , ...,\n",
" 28000.75903421, 24846.06853717, 25603.62271956],\n",
" ...,\n",
" [28532.23134009, 28493.8273213 , 28000.75903421, ...,\n",
" 0. , 14448.29575676, 18710.76846276],\n",
" [25291.03668667, 24995.94592241, 24846.06853717, ...,\n",
" 14448.29575676, 0. , 4264.61589559],\n",
" [25994.5064412 , 25633.14869253, 25603.62271956, ...,\n",
" 18710.76846276, 4264.61589559, 0. ]])\n",
"Coordinates:\n",
" * id (id) <U6 3kB 'ams1' 'ams2' 'ams3' ... 'ams133' 'ams134'\n",
" * id_neighbor (id_neighbor) <U6 3kB 'ams1' 'ams2' ... 'ams133' 'ams134'<xarray.Dataset> Size: 57kB\n",
"Dimensions: (id: 134, n_closest: 25)\n",
"Coordinates:\n",
" * id (id) <U6 3kB 'ams1' 'ams2' 'ams3' ... 'ams133' 'ams134'\n",
"Dimensions without coordinates: n_closest\n",
"Data variables:\n",
" distance (id, n_closest) float64 27kB 0.0 518.8 532.0 ... inf inf inf\n",
" neighbor_id (id, n_closest) object 27kB 'ams1' 'ams2' 'ams3' ... None None"
],
"text/plain": [
"<xarray.DataArray 'neighbor_id' (n_closest: 22)> Size: 176B\n",
"array(['ams63', 'ams61', 'ams70', 'ams58', 'ams71', 'ams74', 'ams51',\n",
" 'ams54', 'ams83', 'ams67', 'ams55', 'ams60', 'ams84', 'ams69',\n",
" 'ams80', 'ams94', 'ams78', 'ams62', 'ams96', 'ams97', 'ams42',\n",
" 'ams41'], dtype=object)\n",
"Coordinates:\n",
" id <U6 24B 'ams63'\n",
"Dimensions without coordinates: n_closest<xarray.Dataset> Size: 167kB\n",
"Dimensions: (cml_id: 345, n_closest: 25)\n",
"Coordinates:\n",
" * cml_id (cml_id) <U21 29kB '10002' '10003' '10004' ... '10363' '10364'\n",
"Dimensions without coordinates: n_closest\n",
"Data variables:\n",
" distance (cml_id, n_closest) float64 69kB 2.745e+03 3.167e+03 ... inf\n",
" neighbor_id (cml_id, n_closest) object 69kB '10001' '10041' ... None None<xarray.DataArray 'neighbor_id' (n_closest: 13)> Size: 104B\n",
"array(['10001', '10041', '10221', '10081', '10101', '10021', '10181',\n",
" '10361', '10121', '10061', '10141', '10341', '10161'], dtype=object)\n",
"Coordinates:\n",
" cml_id <U21 84B '10002'\n",
"Dimensions without coordinates: n_closest