diff --git a/docs/source/howto/create_quicklooks.ipynb b/docs/source/howto/create_quicklooks.ipynb new file mode 100644 index 0000000..98ac46a --- /dev/null +++ b/docs/source/howto/create_quicklooks.ipynb @@ -0,0 +1,267 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1ff5dc6d-caec-4303-8c22-684792c6cfe5", + "metadata": {}, + "source": [ + "# Some dropsonde quicklook examples" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "947b5fe5-e348-4d8c-b6fd-de449fc2aa44", + "metadata": {}, + "outputs": [ + { + "ename": "ModuleNotFoundError", + "evalue": "No module named 'halodrops'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[1], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mhalodrops\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mplotting\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m quicklooks \u001b[38;5;28;01mas\u001b[39;00m ql\n", + "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'halodrops'" + ] + } + ], + "source": [ + "from halodrops.plotting import quicklooks as ql" + ] + }, + { + "cell_type": "markdown", + "id": "6b2d5ef7-5111-4978-99dc-36b6a28fed2c", + "metadata": {}, + "source": [ + "### Test plotting for EUREC4A data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a3d0cd28-5bce-46d9-abe8-5c45cc32c14e", + "metadata": {}, + "outputs": [], + "source": [ + "import eurec4a\n", + "from datetime import datetime, date\n", + "from collections import defaultdict\n", + "from functools import reduce\n", + "\n", + "# Access EUREC4A catalog\n", + "cat = eurec4a.get_intake_catalog(use_ipfs=False)\n", + "\n", + "# Get flight segments\n", + "meta = eurec4a.get_flight_segments()\n", + "\n", + "segments = [{**s,\n", + " \"platform_id\": platform_id,\n", + " \"flight_id\": flight_id\n", + " }\n", + " for platform_id, flights in meta.items()\n", + " for flight_id, flight in flights.items()\n", + " for s in flight[\"segments\"]\n", + " ]\n", + "\n", + "segments_by_segment_id = {s[\"segment_id\"]: s for s in segments}\n", + "segments_ordered_by_start_time = list(sorted(segments, key=lambda s: s[\"start\"]))\n", + "\n", + "circles_Jan24 = [s[\"segment_id\"]\n", + " for s in segments_ordered_by_start_time\n", + " if \"circle\" in s[\"kinds\"]\n", + " and s[\"start\"].date() == date(2020,1,24)\n", + " and s[\"platform_id\"] == \"HALO\"\n", + " ]\n", + "\n", + "first_circle_Jan24 = circles_Jan24[0]\n", + "\n", + "dropsonde_ids = segments_by_segment_id[first_circle_Jan24][\"dropsondes\"][\"GOOD\"]" + ] + }, + { + "cell_type": "markdown", + "id": "bde5f87f-503d-48f2-9cfc-8a4a8c9010a9", + "metadata": {}, + "source": [ + "Load dropsonde dataset and select first circle of January 24, 2020:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5f08a63c-b9e8-4d9e-b48c-ae27983a3920", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "ds = cat.dropsondes.JOANNE.level3.to_dask()\n", + "\n", + "# Select dropsondes from Jan 24 2020\n", + "mask_sondes_first_circle_Jan24 = reduce(lambda a, b: a | b, [ds.sonde_id==d for d in dropsonde_ids])\n", + "mask_sondes_third_circle_Jan24 = reduce(lambda a, b: a | b, [ds.sonde_id==d for d in dropsonde_ids])\n", + "\n", + "ds_sondes_first_circle_Jan24 = ds.isel(sonde_id=mask_sondes_first_circle_Jan24)\n", + "ds_sondes_third_circle_Jan24 = ds.isel(sonde_id=mask_sondes_third_circle_Jan24)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14d71026-054c-410f-8723-b44adcf31f0b", + "metadata": {}, + "outputs": [], + "source": [ + "ds_sondes_first_circle_Jan24" + ] + }, + { + "cell_type": "markdown", + "id": "f01f4463-b18a-4562-ac57-0278a5d27a19", + "metadata": {}, + "source": [ + "Access satellite images:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6de0f717-bd26-4eec-8eba-6f43ee98a57a", + "metadata": {}, + "outputs": [], + "source": [ + "satellite_image = ql.get_satellite_data(ds_flight=ds_sondes_first_circle_Jan24,\n", + " satellite_name=\"goes16\", channel = 13, product=\"ABI-L2-CMIPF\")#,\n", + " #extent = (-62,-48,10,30))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5118c746-9a84-4323-b3bc-c1b3c2c61793", + "metadata": {}, + "outputs": [], + "source": [ + "satellite_image" + ] + }, + { + "cell_type": "markdown", + "id": "b5834ca1-2c9f-47c7-bf7b-0c8ce62eb897", + "metadata": {}, + "source": [ + "## Plot dropsonde launch locations on a map (over a GOES-16 image)\n", + "\n", + "By default the satellite image will be taken at the average launch time from all dropsondes in the plot. By default the colormap is chosen to show the flight altitude, but it can be any variable with dimensions (`sonde_id`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "652d22a6-0ff9-45f8-8b69-8608f18d7ae5", + "metadata": {}, + "outputs": [], + "source": [ + "ql.launch_locations_map(ds_flight=ds_sondes_first_circle_Jan24, save_filepath=\"/home/m/m300931/\", \n", + " satellite_data=satellite_image)\n", + " #extent = (-62,-48,10,30))" + ] + }, + { + "cell_type": "markdown", + "id": "61c0e563-6411-41c1-a21e-d4d23c064537", + "metadata": {}, + "source": [ + "## Plot latitude-time quantities\n", + "\n", + "By default the colormap is chosen to show the flight altitude, but it can be any variable with dimensions (`sonde_id`)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5dc93114-76fc-4ce6-88b9-dab6944428f6", + "metadata": {}, + "outputs": [], + "source": [ + "ql.plot_lat_time(ds_flight=ds_sondes_first_circle_Jan24, save_filepath=\"/home/m/m300931/\")" + ] + }, + { + "cell_type": "markdown", + "id": "e7005a0b-8603-4a8e-8306-af5933ce2a5a", + "metadata": {}, + "source": [ + "## Plot vertical profiles\n", + "The variables and number of plots can be adjusted by the user. By default the quicklook shows temperature, potential temperature, relative humidity, and wind speed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d1315af-7d6e-499e-a1ce-2cda46562724", + "metadata": {}, + "outputs": [], + "source": [ + "ql.plot_profiles(ds_flight=ds_sondes_first_circle_Jan24, save_filepath=\"/home/m/m300931/\")" + ] + }, + { + "cell_type": "markdown", + "id": "c54961a5-873d-4091-b798-d72e42d06a45", + "metadata": {}, + "source": [ + "## Plot dropsonde drift" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "805d161e-72bd-4596-948d-75842b6df1e2", + "metadata": {}, + "outputs": [], + "source": [ + "ql.drift_plots(ds_flight=ds_sondes_first_circle_Jan24, save_filepath=\"/home/m/m300931/\")" + ] + }, + { + "cell_type": "markdown", + "id": "31d1e2b6-0bb0-4b03-b94f-cd0619b7abe7", + "metadata": {}, + "source": [ + "## Saving plots in one pdf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b542d28-3a83-4801-b767-316959ede93d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "HALO-DROPS", + "language": "python", + "name": "my-kernel" + }, + "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.10.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/halodrops/api/plotting.py b/src/halodrops/api/plotting.py new file mode 100644 index 0000000..9c32718 --- /dev/null +++ b/src/halodrops/api/plotting.py @@ -0,0 +1,19 @@ +from halodrops.plotting import quicklooks as ql + + +# Access satellite data +satellite_image = ql.get_satellite_data() + +# Plot launch locations over satellite images +ql.launch_locations_map() + +# Plot longitude/time +ql.plot_lat_time() + +# Plot vertical profiles +ql.plot_profiles() + +# Plot dropsonde drift +ql.drift_plots() + +# Output all quicklooks into a PDF diff --git a/src/halodrops/quicklooks.py b/src/halodrops/quicklooks.py new file mode 100644 index 0000000..852ae36 --- /dev/null +++ b/src/halodrops/quicklooks.py @@ -0,0 +1,369 @@ +from halodrops import sonde +from halodrops.helper import paths +from halodrops.qc import profile + +from importlib import reload +from gogoesgone import processing as pr +from gogoesgone import zarr_access as za + +reload(pr) +reload(za) + +import xarray as xr +import numpy as np +import cartopy as cp +import cartopy.crs as ccrs +from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER +from cartopy.feature import LAND +import cartopy.feature as cfeature +from matplotlib.offsetbox import AnchoredText +import matplotlib.pyplot as plt +import matplotlib.ticker as mticker +from mpl_toolkits.axes_grid1 import make_axes_locatable +import matplotlib.dates as mdates +from datetime import datetime, date +import s3fs +import pandas as pd + + +def convert_time_to_str(time=None, time_format="%Y%m%d %H:%M:%S"): + + """ + Convert input time into desired string format. + """ + + # Ensure time is in correct format + timestamp = (time - np.datetime64("1970-01-01T00:00:00")) / np.timedelta64(1, "s") + datetime_time = datetime.utcfromtimestamp(timestamp) + str_time = datetime_time.strftime(time_format) + + return str_time + + +def get_mean_launch_time(ds_flight=None, time_format="%Y%m%d %H:%M:%S"): + + """ + Compute mean launch time from all sondes in the dataset. + """ + + mean_time = convert_time_to_str(ds_flight.launch_time.mean().values, time_format) + + return mean_time + + +def get_satellite_data( + satellite_time="mean_launch_time", + time_format="%Y%m%d %H:%M:%S", + ds_flight=None, + satellite_name="goes16", + channel=13, + product="ABI-L2-CMIPF", + extent=None, +): + + """ + Access satellite data for nearest time, map to lon/lat grid, and convert to dataset. + By default use the mean launch time from dropsonde dataset. + """ + + # Get correct time for satellite data. Default is mean of all launches + if satellite_time == "mean_launch_time": + use_time = get_mean_launch_time(ds_flight=ds_flight) + else: + use_time = convert_time_to_str(time=satellite_time) + + # Define default extent based on launch locations, or user-defined extent + if extent is not None: + image_lims = extent + else: + image_lims = ( + ds_flight.lon.min() - 1, + ds_flight.lon.max() + 1, + ds_flight.lat.min() - 1, + ds_flight.lat.max() + 1, + ) + + # Get filepath to satellite data at nearest time. + flist = za.nearest_time_url(use_time, time_format, channel, product, satellite_name) + m = za.get_mapper_from_mzz(flist) + + # Select subset of satellite domain + img = pr.Image(m) + subset = img.subset_region_from_latlon_extents(image_lims, unit="degree") + + return subset + + +def launch_locations_map( + ds_flight=None, + satellite_data=None, + save_filepath="/path/to/save/", + color_coding_var="flight_altitude", + color_coding_cmap="magma", + satellite_time=None, + extent=None, + satellite_cmap="Greys", + satellite_vmin=280, + satellite_vmax=300, +): + + """ + Plot dropsonde launch locations, optionally over satellite images. + """ + + fig = plt.figure(figsize=(10, 8)) + + # Define default extent based on launch locations, or user-defined extent + if extent is not None: + image_lims = extent + else: + image_lims = ( + ds_flight.lon.min() - 1, + ds_flight.lon.max() + 1, + ds_flight.lat.min() - 1, + ds_flight.lat.max() + 1, + ) + + ax = plt.axes(projection=ccrs.PlateCarree()) + ax.coastlines(resolution="50m", linewidth=1.5) + ax.set_extent(image_lims, crs=ccrs.PlateCarree()) + + # Plot satellite image + if satellite_data: + sat_im = satellite_data.CMI.isel(t=0).plot( + ax=ax, + x="lon", + y="lat", + cmap=satellite_cmap, + add_colorbar=True, + cbar_kwargs={ + "pad": 0.05, + "extend": "both", + "aspect": 15, + "shrink": 0.7, + "label": f"{satellite_data.CMI.name} / {satellite_data.CMI.units}", + }, + vmin=satellite_vmin, + vmax=satellite_vmax, + zorder=-1, + transform=ccrs.PlateCarree(), + ) + + # Plot flight path + ax.plot( + ds_flight["lon"].isel(alt=-700), + ds_flight["lat"].isel(alt=-700), + c="red", + linestyle=":", + transform=ccrs.PlateCarree(), + zorder=1, + ) + + # Plot launch locations + im_launches = ax.scatter( + ds_flight["lon"].isel(alt=-700), + ds_flight["lat"].isel(alt=-700), + marker="o", + edgecolor="grey", + s=60, + transform=ccrs.PlateCarree(), + c=ds_flight[color_coding_var], + cmap=color_coding_cmap, + ) + + # Assigning axes ticks + xticks = np.arange(-180, 180, 1) + yticks = np.arange(-90, 90, 1) + + # Setting up the gridlines + gl = ax.gridlines( + crs=ccrs.PlateCarree(), + draw_labels=True, + linewidth=1, + color="gray", + alpha=0.2, + linestyle="--", + ) + gl.xlocator = mticker.FixedLocator(xticks) + gl.ylocator = mticker.FixedLocator(yticks) + gl.xformatter = LONGITUDE_FORMATTER + gl.yformatter = LATITUDE_FORMATTER + gl.xlabel_style = {"size": 10, "color": "k"} + gl.ylabel_style = {"size": 10, "color": "k"} + + # Colorbar + g = fig.colorbar( + im_launches, orientation="horizontal", extend="both", aspect=30, pad=0.05 + ) + g.set_label( + f"{ds_flight[color_coding_var].name} / {ds_flight[color_coding_var].units}", + fontsize=12, + ) + plt.tick_params(labelsize=10) + + # Format time stamp + title_time = convert_time_to_str( + time=satellite_data.t[0].values, time_format="%Y-%m-%d %H:%M" + ) + ax.set_title("") + + plt.title( + f"Sondes {ds_flight.sonde_id.values[0]} to {ds_flight.sonde_id.values[-1]} (Satellite Time = {title_time})", + fontsize=12, + pad=10, + ) + + # Save figure + save_filename = f"{save_filepath}launch-locations-{color_coding_var}-{satellite_data.platform_ID}.png" + plt.savefig(save_filename, dpi=300, bbox_inches="tight") + + +def plot_lat_time( + ds_flight, + color_coding_var="flight_altitude", + color_coding_cmap="magma", + save_filepath="/path/to/save/", +): + + """ + Plot spatio-temporal variation (lat v/s time) of selected variable. + """ + + ax = plt.figure(figsize=(12, 4)) + plt.scatter( + ds_flight["launch_time"].values, + ds_flight["lat"].isel(alt=-700).values, + s=90, + c=ds_flight[color_coding_var], + edgecolor="grey", + cmap=color_coding_cmap, + ) + plt.xlim( + np.min(ds_flight["launch_time"].values) - np.timedelta64(4, "m"), + np.max(ds_flight["launch_time"].values) + np.timedelta64(4, "m"), + ) + plt.gca().spines["right"].set_visible(False) + plt.gca().spines["top"].set_visible(False) + g = plt.colorbar() + g.set_label( + f"{ds_flight[color_coding_var].name} / {ds_flight[color_coding_var].units}", + fontsize=12, + ) + + myFmt = mdates.DateFormatter("%H:%M") + plt.gca().xaxis.set_major_formatter(myFmt) + plt.xlabel("Time / UTC", fontsize=12) + plt.ylabel("Latitude / $\degree$N", fontsize=12) + plt.title( + f"Sondes {ds_flight.sonde_id.values[0]} to {ds_flight.sonde_id.values[-1]}", + fontsize=12, + ) + + save_filename = f"{save_filepath}spatiotemporal-variation-{color_coding_var}.png" + + plt.savefig( + save_filename, + dpi=300, + bbox_inches="tight", + ) + + +def plot_profiles( + ds_flight, + r=["ta", "theta", "rh", "wspd", "wdir"], + r_titles=[ + "T / $\degree$C", + "$\\theta$ / K", + "RH / %", + "Wind speed / ms$^{-1}$", + "Wind direction / $\degree$", + ], + row=1, + col=4, + save_filepath="/path/to/save/", +): + + """ + Plot vertical profiles of specified variables. + """ + + f, ax = plt.subplots(row, col, sharey=True, figsize=(12, 6)) + + for j in range(col): + d = ds_flight[r[j]] + for i in range(1, len(ds_flight["launch_time"]) - 1): + ax[j].plot( + d.isel(sonde_id=i), + ds_flight["alt"] / 1000, + c="grey", + alpha=0.25, + linewidth=0.5, + ) + + ax[j].plot( + np.nanmean(d, axis=0), + ds_flight["alt"] / 1000, + linewidth=3, + c="k", + ) + ax[j].set_xlabel(r_titles[j], fontsize=12) + ax[j].spines["right"].set_visible(False) + ax[j].spines["top"].set_visible(False) + if j == 0: + ax[j].set_ylabel("Altitude / km", fontsize=12) + + plt.suptitle( + f"Sondes {ds_flight.sonde_id.values[0]} to {ds_flight.sonde_id.values[-1]}", + fontsize=12, + ) + + save_filename = f"{save_filepath}vertical-profiles-measured-quantities.png" + + plt.savefig(save_filename, dpi=300, bbox_inches="tight") + + +def drift_plots(ds_flight=None, save_filepath="/path/to/save/"): + + print("Plotting drift in lat and lon...") + + f, ax = plt.subplots(1, 2, sharey=True, figsize=(10, 5)) + + for i in range(len(ds_flight["launch_time"])): + + max_id = np.max(np.where(~np.isnan(ds_flight["lon"].isel(sonde_id=i)))) + + ax[0].plot( + ds_flight["lat"].isel(sonde_id=i) + - ds_flight["lat"].isel(sonde_id=i).isel(alt=max_id), + ds_flight["alt"] / 1000, + linewidth=1.5, + c="grey", + alpha=0.75, + ) + + ax[0].set_xlabel("Drift in Latitude / $\degree$", fontsize=12) + ax[0].set_ylabel("Altitude / km", fontsize=12) + ax[0].spines["right"].set_visible(False) + ax[0].spines["top"].set_visible(False) + + ax[1].plot( + ds_flight["lon"].isel(sonde_id=i) + - ds_flight["lon"].isel(sonde_id=i).isel(alt=max_id), + ds_flight["alt"] / 1000, + linewidth=1.5, + c="grey", + alpha=0.75, + ) + + ax[1].set_xlabel("Drift in Longitude / $\degree$", fontsize=12) + ax[1].spines["right"].set_visible(False) + ax[1].spines["top"].set_visible(False) + + plt.suptitle( + f"Sondes {ds_flight.sonde_id.values[0]} to {ds_flight.sonde_id.values[-1]}", + fontsize=12, + ) + + save_filename = f"{save_filepath}drift-in-lat-lon.png" + + plt.savefig(save_filename, dpi=300, bbox_inches="tight")