Trixels truncated at dateline #143
Description
Here's an example program to demonstrate this.
#! /usr/bin/env python -tt
# -*- coding: utf-8; mode: python -*-
r"""
dateline_issue.py
~~~~~~~~~~~~~~~~~
"""
# Standard Imports
import os
import pickle
# Third-Party Imports
from matplotlib import cm
from matplotlib import font_manager
from shapely.geometry import Polygon
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import dask
import geopandas
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas
# STARE Imports
import pystare
import starepandas
# ------------------------------------------------------------------------------
def main():
##
# Tweak Pandas screen dump
pandas.set_option("display.max_columns", None)
pandas.set_option("display.width", 500)
pandas.set_option("display.colheader_justify", "right")
pandas.set_option("display.precision", 3)
pandas.options.mode.copy_on_write = True
##
# Tweak GeoPandas screen dump
geopandas.options.display_precision = 3
geopandas.options.use_pygeos = False
# Mike's local copy
# feature_data = "/Users/mbauer/tmp/data/tablespace/xcal/pickles/featuredb.pkl"
# imerg_sidecar = "IMERG_stare.nc"
# for some reason the FlexFS servers /bayesics folder is empty, but was something like this
feature_data = "/tablespace/xcal/pickles/featuredb.pickle"
imerg_sidecar = "/bayesics/p4/Source/GPM/IMERG_stare.nc"
draw_plots = [False, True][0]
##
# Read the IMERG Features file into a STAREPandas DataFrame
with open(feature_data, 'rb') as f:
imerg_sdf = pickle.load(f)
##
# Extract a single feature/event
# subset_sdf = <class 'starepandas.staredataframe.STAREDataFrame'>
"""
index label timestamp x y cell_areas tot_area precips tot_precip sids cover trixels
0 6077 37 2021-02-05 14:30:00 [542, 542, 542, 543, 543, 544, 544] [3233, 3234, 3235, 3234, 3235, 3233, 3234] [100345481.30079782, 100345481.30079782, 10034... 7.032e+08 [1.355738, 2.0544746, 1.5360532, 1.3778167, 1.... 5.003e+05 [904835998274304745, 904828974740792201, 90481... [904812307752681481, 904814506775937033, 90481... MULTIPOLYGON (((143.676 35.644, 143.460 35.663...
1 6078 37 2021-02-05 15:00:00 [541, 541, 542, 542, 543, 543, 543, 544] [3234, 3235, 3234, 3235, 3233, 3234, 3235, 3233] [100219248.63592072, 100219248.63592072, 10034... 8.031e+08 [1.2717158, 1.9754196, 1.1044679, 1.495765, 1.... 5.529e+05 [904827950646811209, 904933268391128745, 90482... [904812307752681481, 904816705799192585, 90482... MULTIPOLYGON (((143.676 35.644, 143.460 35.663...
2 6079 37 2021-02-05 15:30:00 [540, 540, 541, 541, 542, 542, 542, 543] [3235, 3236, 3235, 3236, 3234, 3235, 3236, 3234] [100092710.68594712, 100092710.68594712, 10021... 8.021e+08 [1.5025604, 2.3538067, 1.3093063, 1.7862825, 1... 6.881e+05 [904919807797595113, 904934765162898217, 90493... [904812307752681481, 904827700915470345, 90483... MULTIPOLYGON (((143.676 35.644, 143.460 35.663...
3 6080 37 2021-02-05 16:00:00 [534, 534, 534, 539, 539, 539, 539, 540, 540, ... [3230, 3231, 3232, 3236, 3237, 3238, 3239, 323... [99327093.64970006, 99327093.64970006, 9932709... 1.499e+09 [2.548102, 3.5231876, 1.6926019, 1.8397847, 1.... 1.545e+06 [904296384236998153, 904297886683378313, 90429... [904293338264371209, 904295537287626761, 90429... MULTIPOLYGON (((143.317 36.440, 143.172 36.643...
4 6081 37 2021-02-05 16:30:00 [532, 533, 533, 533, 534, 538, 538, 538, 538, ... [3231, 3231, 3232, 3233, 3232, 3238, 3239, 324... [99069463.27244306, 99198429.54881625, 9919842... 1.796e+09 [1.42, 2.77, 3.83, 1.8399999, 1.0699999, 2.0, ... 1.877e+06 [904308668621248745, 904300510651838537, 90430... [904293338264371209, 904299935334137865, 90430... MULTIPOLYGON (((143.317 36.440, 143.172 36.643...
.. ... ... ... ... ... ... ... ... ... ... ... ...
199 6276 37 2021-02-09 18:00:00 [568, 568, 568, 569, 569, 569, 569, 569, 570, ... [3597, 3598, 3599, 3595, 3596, 3597, 3598, 359... [103519134.98180684, 103519134.98180684, 10351... 1.392e+10 [1.1695099, 1.6371919, 1.568149, 1.2314502, 1.... 1.119e+07 [867223753433970025, 867209926233394793, 86714... [867083665757175816, 867092461850198024, 86710... MULTIPOLYGON (((179.472 32.661, 179.834 32.486...
200 6277 37 2021-02-09 18:30:00 [571, 572, 572, 573, 573, 573, 573, 574, 574, ... [3599, 3598, 3599, 3595, 3597, 3598, 3599, 359... [103871730.60377292, 103988630.13166831, 10398... 1.282e+10 [1.9399999, 1.0699999, 2.1299999, 1.0, 1.59999... 1.045e+07 [867208080067029545, 867101236360906281, 86709... [867083665757175816, 867092461850198025, 86709... MULTIPOLYGON (((179.472 32.661, 179.834 32.486...
201 6278 37 2021-02-09 19:00:00 [573, 574, 574, 575, 575, 575, 576, 576, 576, ... [3599, 3597, 3599, 3597, 3598, 3599, 3597, 359... [104105212.892217, 104221478.53028575, 1042214... 1.020e+10 [1.4010634, 1.1362529, 1.8140539, 1.2430454, 1... 1.036e+07 [868931415636743593, 867116255144589001, 86893... [867085864780431369, 867114452082753545, 86778... MULTIPOLYGON (((179.727 32.334, 179.653 32.574...
202 6279 37 2021-02-09 19:30:00 [576, 577, 577, 578, 578, 579, 579, 579, 580, ... [3599, 3598, 3599, 3598, 3599, 3597, 3598, 359... [104453057.02328122, 104568369.17277576, 10456... 3.893e+09 [1.2154739, 1.3653224, 1.4656936, 1.0716105, 1... 3.633e+06 [868943611143033993, 868920339897876713, 86894... [867787353198952456, 867796149291974665, 86780... MULTIPOLYGON (((179.981 32.006, 179.769 31.702...
203 6280 37 2021-02-09 20:00:00 [584, 584, 585, 585, 586] [3598, 3599, 3598, 3599, 3599] [105366615.73733008, 105366615.73733008, 10547... 5.273e+08 [1.0232406, 1.1071403, 1.07733, 1.2315224, 1.2... 2.991e+05 [867801516070089001, 867803219036514633, 86779... [867796149291974665, 867800547338485769, 86780... MULTIPOLYGON (((179.915 31.223, 179.842 31.463...
[204 rows x 12 columns]
"""
subset_sdf = imerg_sdf[imerg_sdf["label"].isin([37])]
subset_sdf = subset_sdf.reset_index()
del imerg_sdf
# ##
# # TMP save time debugging (comment out to "Find the SIDs for verts/geometry" after save)
# with open('tmp.pkl', "wb") as pfile:
# pickle.dump(subset_sdf, pfile, pickle.HIGHEST_PROTOCOL)
# ##
# # TMP save time debugging (comment out above to "Find the SIDs for verts/geometry")
# with open('tmp.pkl', 'rb') as f:
# subset_sdf = pickle.load(f)
# <class 'pandas.core.series.Series'>
single_sdf = subset_sdf.iloc[126, :]
print(f"{type(single_sdf) = }")
"""
sids (4060) : [874279450188275529 ... 3584450484248280905]
cover (617) : [874278869849341960 ... 3584449687991615497]
tixels : MULTIPOLYGON (with 618 Polygons)
Looping over these trixel Polygons I see that no trixels have longitudes past -180.
tmp[ 0] max/min lon = +168.44 ... +168.84
tmp[ 1] max/min lon = +168.74 ... +168.94
tmp[ 2] max/min lon = +168.74 ... +168.94
...
tmp[615] max/min lon = +176.25 ... +176.59
tmp[616] max/min lon = +176.00 ... +176.34
tmp[617] max/min lon = -180.00 ... -179.97
Min Lon -180.00
Max Lon -179.97
Now if I convert the SIDs into vertices I see something different:
verts = starepandas.to_vertices(sids, wrap_lon=False)
corner lats (12180, 30.74, 48.46):
corner lons (12180, 162.82, 180.03):
center lats (4060, 30.87, 48.37):
center lons (4060, 162.93, 179.89):
So the corner lons just barely cross the DL at 180.03 deg.
But if I switch wrap_lon:
verts = starepandas.to_vertices(sids, wrap_lon=True)
corner lats (12180, 30.74, 48.46):
corner lons (12180, -179.97, 179.92):
center lats (4060, 30.87, 48.37):
center lons (4060, 162.93, 179.89):
Now the corners lons span the interior between the DL +/-179 deg.
"""
sids = single_sdf["sids"]
cover = single_sdf["cover"]
tixels = single_sdf["trixels"]
print(f"sids {sids.shape} : [{np.amin(sids)} ... {np.amax(sids)}]")
print(f"cover {cover.shape}: [{np.amin(cover)} ... {np.amax(cover)}]")
print(f"tixels : MULTIPOLYGON(with {len(tixels.geoms)} Polygons)")
lons = []
for inner_poly in range(len(tixels.geoms)):
tmp = list(tixels.geoms[inner_poly].exterior.coords)
lons.extend([_[0] for _ in tmp])
tmp_lons = sorted([_[0] for _ in tmp])
print(f"tmp[{inner_poly:3d}] max/min lon = {tmp_lons[0]:+8.2f} ... {tmp_lons[-1]:+8.2f}")
print(f"Min lon {min(tmp_lons):+8.2f}")
print(f"Max lon {max(tmp_lons):+8.2f}")
verts = starepandas.to_vertices(sids, wrap_lon=False)
print(f"\n{verts = }")
plabs = ("corner lats", "corner lons", "center lats", "center lons")
for pidx, part in enumerate(verts):
print(f"\t{plabs[pidx]} ({len(part)}, {np.amin(part):10.2f}, {np.amax(part):10.2f}): {part}")
verts = starepandas.to_vertices(sids, wrap_lon=True)
print(f"\n{verts = }")
plabs = ("corner lats", "corner lons", "center lats", "center lons")
for pidx, part in enumerate(verts):
print(f"\t{plabs[pidx]} ({len(part)}, {np.amin(part):10.2f}, {np.amax(part):10.2f}): {part}")
if draw_plots:
plot_dpi = 300
plot_size = [12, 4]
# WGS-84 Earth equatorial radius at sea level (meters)
# STARE uses ccrs.Globe(datum='WGS84', ellipse='WGS84') https://epsg.io/4326
globe = ccrs.Globe(datum='WGS84', ellipse='WGS84')
lon_0_global = 0.0
lat_0_global = 90.0
##
# Map extent
min_lat, max_lat = 20, 70
min_lon, max_lon = -180, 180
# Determine the map extent (left_lon, right_lon, lower_lat, upper_lat) for the plot domain
map_extent = (min_lon, max_lon, min_lat, max_lat)
##
# Utility CRS
# Geodetic:
# A 3D/spherical CRS based on latitude and longitude where geographical distance and coordinates are measured in degrees.
geod_crs = ccrs.Geodetic(globe=globe)
# PlateCarree Projection (AKA Equirectangular or Equidistant Cylindrical):
# A 2D CRS with a flat topology and Euclidean distance (m).
flat_crs = ccrs.PlateCarree(central_longitude=lon_0_global, globe=globe)
##
# Main Map CRS
# PlateCarree Projection (AKA Equirectangular or Equidistant Cylindrical):
# A 2D CRS with a flat topology and Euclidean distance (m).
map_crs = ccrs.PlateCarree(central_longitude=lon_0_global, globe=globe)
proj_title = "PlateCarree"
# Lambert Azimuthal Equal-Area
# A 2D CRS centered on (map_lat0, map_lon0), which is then the zero point of the flat topology w/ Euclidean distance (m).
# polar_crs = ccrs.LambertAzimuthalEqualArea(central_longitude=lon_0_global, central_latitude=lat_0_global, globe=globe)
# Tweak CRS
map_crs._threshold /= 100.
for index, row in subset_sdf.iterrows():
#if index < 126:
if index < 146:
continue
##
# Make Figure
# --------------------------------------------------------------------------
fig = plt.figure(figsize=(12, 4), frameon=True)
geo_axes = plt.axes(projection=map_crs)
##
# Set map extent
geo_axes.set_xlim(left=min_lon, right=max_lon)
geo_axes.set_ylim(bottom=min_lat, top=max_lat)
##
# Color Land/Sea
geo_axes.add_feature(cfeature.LAND, alpha=0.5)
geo_axes.add_feature(cfeature.COASTLINE, alpha=0.5)
dlat = 10
par = np.arange(min_lat, max_lat + dlat, dlat)
dlon = 30
mer = np.arange(min_lon, max_lon + dlon, dlon)
gl = geo_axes.gridlines(crs=map_crs, draw_labels=True, linewidth=0.5, color='k', alpha=0.5,
linestyle='--', x_inline=False, y_inline=False)
sdf_ = starepandas.STAREDataFrame({'index': [row.index]}, sids=[row.sids], trixels=[row.trixels])
sdf_.plot(ax=geo_axes, trixels=True, boundary=True, figsize=plot_size, aspect=None, zorder=2, linewidth=0.5, color='r', alpha=0.7, transform=geod_crs)
geo_axes.set_title(f"Feature {str(row.label)}: index {str(index)}")
plt.show()
plt.clf()
plt.close('all')
return
# ---Start of main code block.
if __name__ == '__main__':
#
# Local DASK Processes (https://dask.pydata.org/en/latest/scheduling.html?highlight=multiprocessing)
#
# To use DASK locally (i.e., cores on the machine you are working on) your python script must have a main()
# method and you need to import dask and issue the following directive. If not, STARE requests such as
# sids = starepandas.sids_from_geoseries(polys, level=10, convex=False, force_ccw=True, n_partitions=2, num_workers=4)
# Throws a RuntimeError:
# An attempt has been made to start a new process before the
# current process has finished its bootstrapping phase.
#
# This probably means that you are not using fork to start your
# child processes and you have forgotten to use the proper idiom
# in the main module:
#
# if __name__ == '__main__':
# freeze_support()
# ...
#
# The "freeze_support()" line can be omitted if the program
# is not going to be frozen to produce an executable.
#
dask.config.set(scheduler='processes')
##
# Run the main routine
main()I should add some potential dependency issues that might be at work here (e.g., Pandas 2.0).
| Name | Version |
|---|---|
| astropy | 5.2.2 |
| cartopy | 0.21.1 |
| gdal | 3.6.3 |
| geopandas | 0.12.2 |
| geopandas-base | 0.12.2 |
| geos | 3.11.2 |
| h5py | 3.8.0 |
| hdf4 | 4.2.15 |
| hdf5 | 1.12.2 |
| hdfeos2 | 2.20 |
| matplotlib | 3.7.1 |
| numpy | 1.24.2 |
| pandas | 2.0.0 |
| proj | 9.1.1 |
| pygeos | 0.14 |
| pyhdf | 0.10.5 |
| pyproj | 3.5.0 |
| pyshp | 2.3.1 |
| pytest | 7.3.1 |
| python | 3.11.3 |
| shapely | 2.0.1 |
| STARE | Installs |
| pystare | 0.8.12 STARE |
| staremaster | 0.0.4 STARE |
| starepandas | 0.6.6 STARE |