Sample satellite data

Some graphical tests of current primary functions using the sample satellite data.

import warnings

import cartopy
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import tams

warnings.filterwarnings("ignore", category=cartopy.io.DownloadWarning)

%matplotlib inline

xr.set_options(display_expand_data=False)

Show code cell output Hide code cell output

ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/tams/conda/latest/share/proj failed

<xarray.core.options.set_options at 0x7fe6cd1b9fc0>

Load data

tb = tams.load_example_tb().isel(time=slice(4))
tb

<xarray.DataArray 'ch9' (time: 4, y: 714, x: 2829)>
295.9 295.9 295.8 295.8 295.8 295.9 295.9 295.8 ... nan nan nan nan nan nan nan
Coordinates:
    lon      (time, y, x) float32 ...
    lat      (time, y, x) float32 ...
  * time     (time) datetime64[ns] 2006-09-01 ... 2006-09-01T06:00:00
Dimensions without coordinates: y, x
Attributes:
    units:      K
    long_name:  Brightness temperature

xarray.DataArray

'ch9'

• time: 4
• y: 714
• x: 2829

• 295.9 295.9 295.8 295.8 295.8 295.9 295.9 ... nan nan nan nan nan nan

  array([[[295.92957, 295.92957, 295.8033 , ..., 292.85956, 292.85956,
           292.7298 ],
          [295.8033 , 295.8033 , 295.92957, ..., 293.11865, 292.98917,
           292.46976],
          [295.8033 , 295.67694, 295.8033 , ..., 292.85956, 292.98917,
           292.59985],
          ...,
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan]],
  
         [[295.8033 , 295.8033 , 295.8033 , ..., 287.54327, 285.89505,
           284.92136],
          [296.05566, 295.8033 , 295.4237 , ..., 288.76266, 287.67944,
           286.0334 ],
          [295.92957, 295.8033 , 295.55038, ..., 290.1013 , 289.70148,
           288.08698],
  ...
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan]],
  
         [[294.14902, 293.50616, 294.14902, ..., 290.7644 , 289.03174,
           289.56784],
          [294.53293, 293.8923 , 294.14902, ..., 290.8965 , 289.43405,
           290.1013 ],
          [294.53293, 294.27713, 294.4051 , ..., 290.8965 , 290.36703,
           291.02853],
          ...,
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan],
          [      nan,       nan,       nan, ...,       nan,       nan,
                 nan]]], dtype=float32)

• Coordinates: (3)
  • lon
    (time, y, x)
    float32
    ...

    _CoordinateAxisType :

    Lon

    standard_name :

    longitude

    valid_max :

    180.0

    long_name :

    Longitude

    valid_min :

    -180.0

    units :

    degrees_east

    [8079624 values with dtype=float32]

  • lat
    (time, y, x)
    float32
    ...

    _CoordinateAxisType :

    Lat

    standard_name :

    latitude

    valid_max :

    90.0

    long_name :

    Latitude

    valid_min :

    -90.0

    units :

    degrees_north

    [8079624 values with dtype=float32]

  • time
    (time)
    datetime64[ns]
    2006-09-01 ... 2006-09-01T06:00:00

    array(['2006-09-01T00:00:00.000000000', '2006-09-01T02:00:00.000000000',
           '2006-09-01T04:00:00.000000000', '2006-09-01T06:00:00.000000000'],
          dtype='datetime64[ns]')

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2006-09-01 00:00:00', '2006-09-01 02:00:00',
                   '2006-09-01 04:00:00', '2006-09-01 06:00:00'],
                  dtype='datetime64[ns]', name='time', freq='2H'))

• Attributes: (2)

  units :

  K

  long_name :

  Brightness temperature

tb.isel(time=0).plot(x="lon", y="lat", size=5, aspect=2.5);

Identify cloud elements (CEs)

times = tb.time
contour_sets, contour_sets_219 = tams.identify(tb)
contour_sets[0].head()

	geometry	area_km2	inds219	area219_km2	cs219
0	POLYGON ((43.10657 10.95914, 43.15180 10.96099...	75505.811057	[2, 3, 4, 5, 6]	29765.307905	MULTIPOLYGON (((43.11648 11.39076, 43.32919 11...
1	POLYGON ((35.04399 10.35363, 35.05942 10.35526...	180394.270027	[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]	56096.842040	MULTIPOLYGON (((35.26395 11.24077, 35.30297 11...
2	POLYGON ((26.33872 0.00000, 26.89784 0.00000, ...	17858.344651	[22]	8989.069746	MULTIPOLYGON (((26.59017 0.02994, 26.78338 0.0...
3	POLYGON ((22.11254 3.49762, 23.49709 3.55631, ...	99926.168432	[25, 26]	76675.055415	MULTIPOLYGON (((22.08261 3.51624, 22.11310 3.5...
4	POLYGON ((12.96943 2.84170, 13.53301 2.84443, ...	378232.478571	[31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 4...	109994.951561	MULTIPOLYGON (((16.02669 5.74587, 16.04421 5.7...

contour_sets_219[0].head()

	geometry	area_km2
0	POLYGON ((46.45179 11.43426, 46.90813 11.48031...	2861.636356
1	POLYGON ((46.41562 12.01331, 46.46724 12.02231...	1005.761023
2	POLYGON ((43.11648 11.39076, 43.32919 11.52867...	376.874022
3	POLYGON ((43.18207 11.28760, 43.22879 11.29522...	21270.230923
4	POLYGON ((41.55725 11.69260, 41.57923 11.69670...	44.893865

# Simple plot to test 219 matching
m, n = 0, 1  # time, contour #
fig, ax = plt.subplots()
c = contour_sets[m].iloc[[n]]
c.plot(ax=ax)
c.cs219.plot(color="red", ax=ax, alpha=0.4);

Track CE groups between times

cs = tams.track(contour_sets, times, u_projection=-5).reset_index(drop=True)
cs.head()

	geometry	area_km2	inds219	area219_km2	cs219	time	itime	dtime	mcs_id
0	POLYGON ((43.10657 10.95914, 43.15180 10.96099...	75505.811057	[2, 3, 4, 5, 6]	29765.307905	MULTIPOLYGON (((43.11648 11.39076, 43.32919 11...	2006-09-01	0	0 days 02:00:00	0
1	POLYGON ((35.04399 10.35363, 35.05942 10.35526...	180394.270027	[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]	56096.842040	MULTIPOLYGON (((35.26395 11.24077, 35.30297 11...	2006-09-01	0	0 days 02:00:00	1
2	POLYGON ((26.33872 0.00000, 26.89784 0.00000, ...	17858.344651	[22]	8989.069746	MULTIPOLYGON (((26.59017 0.02994, 26.78338 0.0...	2006-09-01	0	0 days 02:00:00	2
3	POLYGON ((22.11254 3.49762, 23.49709 3.55631, ...	99926.168432	[25, 26]	76675.055415	MULTIPOLYGON (((22.08261 3.51624, 22.11310 3.5...	2006-09-01	0	0 days 02:00:00	3
4	POLYGON ((12.96943 2.84170, 13.53301 2.84443, ...	378232.478571	[31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 4...	109994.951561	MULTIPOLYGON (((16.02669 5.74587, 16.04421 5.7...	2006-09-01	0	0 days 02:00:00	4

tams.plot_tracked(cs)

# CTT movie -- WIP

import cartopy.crs as ccrs
import regionmask

i = -1  # select time
size = 3
cx = cy = 3

tb_ = tb.isel(time=i)
if cx > 1 or cy > 1:
    tb_ = tb_.coarsen(x=cx, y=cy, boundary="trim").mean()

x0, x1 = tb_.lon.values.min(), tb_.lon.values.max()
y0, y1 = tb_.lat.values.min(), tb_.lat.values.max()

aspect = (x1 - x0) / (y1 - y0)
proj = ccrs.Mercator()
tran = ccrs.PlateCarree()

fig = plt.figure(figsize=(size * aspect, size + 1))
gs = fig.add_gridspec(
    2, 2,
    width_ratios=(1, 1), height_ratios=(aspect * 2 + 1, 1),
    left=0.1, right=0.9, bottom=0.1, top=0.9,
    wspace=0.05, hspace=0.18,
)

ax = fig.add_subplot(gs[0, :], projection=proj)
ax.set_extent([x0, x1, y0, y1], crs=tran)
ax.gridlines(draw_labels=True)

ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[1, 1])

# Background -- CTT
tb_.plot(
    x="lon", y="lat",
    cmap="gray_r", ax=ax, cbar_ax=ax2,
    transform=tran,
    cbar_kwargs=dict(orientation="horizontal"),
)

# CEs with colored precip
shapes = cs[["geometry"]]
regions = regionmask.from_geopandas(shapes)
print("masking")
mask = regions.mask(tb_)
masked = tb_.where(mask >= 0)
print("plotting masked")
masked.plot.pcolormesh(
    x="lon", y="lat",
    ax=ax, cbar_ax=ax3, transform=tran, alpha=0.6,
    cbar_kwargs=dict(orientation="horizontal"),
)

# Tracks up to this time
print("plotting tracks")
t = tb_.time.values
assert np.isscalar(t)
for _, g in cs.groupby("mcs_id"):
    g_ = g[g.time <= t]
    c = g_.to_crs("EPSG:32663").centroid.to_crs("EPSG:4326")
    ax.plot(c.x, c.y, ".-", c="r", lw=2, transform=tran)

masking

plotting masked
plotting tracks

Classify

cs = tams.classify(cs)
cs.head()

	geometry	area_km2	inds219	area219_km2	cs219	time	itime	dtime	mcs_id	mcs_class
0	POLYGON ((43.10657 10.95914, 43.15180 10.96099...	75505.811057	[2, 3, 4, 5, 6]	29765.307905	MULTIPOLYGON (((43.11648 11.39076, 43.32919 11...	2006-09-01	0	0 days 02:00:00	0	DSL
1	POLYGON ((35.04399 10.35363, 35.05942 10.35526...	180394.270027	[10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]	56096.842040	MULTIPOLYGON (((35.26395 11.24077, 35.30297 11...	2006-09-01	0	0 days 02:00:00	1	MCC
2	POLYGON ((26.33872 0.00000, 26.89784 0.00000, ...	17858.344651	[22]	8989.069746	MULTIPOLYGON (((26.59017 0.02994, 26.78338 0.0...	2006-09-01	0	0 days 02:00:00	2	DSL
3	POLYGON ((22.11254 3.49762, 23.49709 3.55631, ...	99926.168432	[25, 26]	76675.055415	MULTIPOLYGON (((22.08261 3.51624, 22.11310 3.5...	2006-09-01	0	0 days 02:00:00	3	MCC
4	POLYGON ((12.96943 2.84170, 13.53301 2.84443, ...	378232.478571	[31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 4...	109994.951561	MULTIPOLYGON (((16.02669 5.74587, 16.04421 5.7...	2006-09-01	0	0 days 02:00:00	4	DLL