diff --git a/notebooks/data_subsets.ipynb b/notebooks/data_subsets.ipynb new file mode 100644 index 0000000..79c2872 --- /dev/null +++ b/notebooks/data_subsets.ipynb @@ -0,0 +1,215 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "aa472915-d85c-4358-b686-48eb52f4075f", + "metadata": {}, + "outputs": [], + "source": [ + "# import packages:\n", + "import numpy as np # Numeric Python\n", + "import matplotlib.pyplot as plt # Plotting routines\n", + "import h5py # general HDF5 reading/writing library\n", + "import rioxarray as rx # Package to read raster data from hdf5 files\n", + "from pyproj import Transformer, CRS # libraries to allow coordinate transforms\n", + "import glob # Package to locate files on disk\n", + "import os # File-level utilities\n", + "import re # regular expressions for string interpretation\n", + "import icepyx as ipx # Package to interact with ICESat-2 online resources\n", + "from sliderule import icesat2 # Package for online ICESat-2 processing" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "581f3be5-a81e-4455-b7c3-66086a5a428f", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd59afa5-d9dc-4ee8-85f1-7a9e634b1b24", + "metadata": {}, + "outputs": [], + "source": [ + "# logins, etc.\n", + "\n", + "#HOST = 'https://urs.earthdata.nasa.gov'\n", + "#ipx.core.Earthdata.Earthdata('ben_smith','whatever@whatever.io', HOST).login()\n", + "\n", + "url=\"icesat2sliderule.org\"\n", + "icesat2.init(url, verbose=False)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ba30c90-37de-43c4-9f88-16fc412ef500", + "metadata": {}, + "outputs": [], + "source": [ + "# Annika's bounding box:\n", + "# x,y\n", + "#-340,-80\n", + "#-480,-170\n", + "\n", + "XR= np.array([-480, -340])*1.e3\n", + "YR= np.array([-170, -80])*1.e3\n", + "\n", + "# shrink down to a tiny box in the center:\n", + "XR=np.mean(XR)+np.array([-5.e3, 5.e3])\n", + "YR=np.mean(YR)+np.array([-5.e3, 5.e3])\n", + "\n", + "\n", + "# Prepare coordinate transformations between lat/lon and polar stereographic\n", + "crs=CRS.from_epsg(3031)\n", + "to_xy_crs=Transformer.from_crs(crs.geodetic_crs, crs)\n", + "to_geo_crs=Transformer.from_crs(crs, crs.geodetic_crs)\n", + "\n", + "corners_lat, corners_lon=to_geo_crs.transform(np.array(XR)[[0, 1, 1, 0, 0]], np.array(YR)[[0, 0, 1, 1, 0]])\n", + "latlims=[np.min(corners_lat), np.max(corners_lat)]\n", + "lonlims=[np.min(corners_lon), np.max(corners_lon)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d183a6e7-b26a-4403-8026-17f09ac8cbcc", + "metadata": {}, + "outputs": [], + "source": [ + "# run a slideRule ATL06 query. Just ask for cycle 8 (Antarctic winter, 2020)\n", + "# to avoid getting swamped right away\n", + "\n", + "# See parameters here:\n", + "# http://icesat2sliderule.org/rtd/user_guide/ICESat-2.html\n", + "params= { 'poly':[{'lon':this_lon, 'lat':this_lat} for this_lon, this_lat in zip(corners_lon, corners_lat)],\n", + " 'srt':3,\n", + " 'cnf':1,\n", + " 'len':10,\n", + " 'res':10,\n", + " 'ats':5,\n", + " 'cnt':10,\n", + " 'cycle':8,\n", + " 'maxi': 10,\n", + " 'pass_invalid':False}\n", + "\n", + "D_IS_SR=icesat2.atl06p(params, \n", + " asset=\"nsidc-s3\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e8aa520b-a61b-44f8-a9ac-214ba87a5d66", + "metadata": {}, + "outputs": [], + "source": [ + "D_IS_SR[0:3]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "40262c5d-9842-494e-8f94-9a95f722cf59", + "metadata": {}, + "outputs": [], + "source": [ + "lon=np.array([gi.x for gi in D_IS_SR.geometry])\n", + "lat=np.array([gi.y for gi in D_IS_SR.geometry])\n", + "xy=to_xy_crs.transform(lat, lon)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6eaf0a5f-57f0-4565-acb1-ffa1e1f01dd8", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(); plt.scatter(xy[0], xy[1], c=D_IS_SR['h_mean'], cmap='Spectral'); plt.colorbar(label='surface height')\n", + "\n", + "\n", + "plt.figure(); plt.scatter(xy[0], xy[1], c=D_IS_SR['rms_misfit'], vmin=0.1, vmax=1); plt.colorbar(label='rms_misfit')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27a8b80c-aa9a-4639-8e44-e9ebc8bcc2c8", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(); plt.hist(D_IS_SR['rms_misfit'], np.arange(0, 2, 0.01));" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef6c1b6e-2dc7-422a-9be4-0c9eb46c701f", + "metadata": {}, + "outputs": [], + "source": [ + "sigma_extra=np.sqrt(np.maximum(0, D_IS_SR['rms_misfit']**2-(0.68e-9*1.5e8)**2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf81d091-8aaf-49b3-83c9-8b38ad04b0a1", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(); plt.hist(sigma_extra, np.arange(0, .5, 0.01));" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a58b03a6-1f54-4a67-8848-870bf152fbff", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(); \n", + "\n", + "plt.scatter(xy[0], xy[1], c=sigma_extra, vmin=0, vmax=0.5); plt.colorbar(label='roughness?')\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e16f3565-d366-4fcb-9389-298bfdab2257", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}