<jupyter_start><jupyter_text>Eddies identification in oceanographic stations using geospatial analisys*****Author:** [<NAME>](https://github.com/andresawa/)  Laboratório de Estudos dos Oceanos e Clima – LEOC, Instituto de Oceanografia, Universidade Federal do Rio Grande.  email: <EMAIL>  **Last change:** 12/10/2021***  **Part 1: Processing** Import libraries<jupyter_code>import numpy as np
import pandas as pd
import xarray as xr
import geopandas as gpd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from shapely.geometry import Point<jupyter_output><empty_output><jupyter_text>Define the temporal and geographic limits<jupyter_code>lon_min, lon_max, lat_min, lat_max, time_min, time_max = -40, 17, -31, -19, np.datetime64("2009-10-24"), np.datetime64("2009-12-20")<jupyter_output><empty_output><jupyter_text>Load the eddies database (META3.1exp) and subset by position/time<jupyter_code># Cyclonic eddies
# Load the eddies and subset using the limits (latitude and time)
cyc = xr.open_dataset("META3_1_cyclonic.nc")
cyc_subset = cyc.sel(obs = 
                     (cyc.latitude >= lat_min) & 
                     (cyc.latitude <= lat_max) & 
                     (cyc.time >= time_min) & 
                     (cyc.time <= time_max)
                    )

# Transform the longitudes > 180 to negative values (change coordinates from [0 ~ 360] to [-180 ~ 180])
cyc_subset["longitude"] = (cyc_subset.longitude + 180) % 360 - 180

# Filter eddies by longitude
cyc_filter = cyc_subset.sel(obs = 
                            (cyc_subset.longitude >= lon_min) & 
                            (cyc_subset.longitude <= lon_max)  
                           )
# Anticyclonic eddies
# Load the eddies and subset using the limits (latitude and time)
acyc = xr.open_dataset("META3_1_anticyclonic.nc")
acyc_subset = acyc.sel(obs = 
                       (acyc.latitude >= lat_min) & 
                       (acyc.latitude <= lat_max) & 
                       (acyc.time >= time_min) & 
                       (acyc.time <= time_max)
                      )

# Transform the longitudes > 180 to negative values (change coordinates from [0 ~ 360] to [-180 ~ 180])
acyc_subset["longitude"] = (acyc_subset.longitude + 180) % 360 - 180

# Filter eddies by longitude
acyc_filter = acyc_subset.sel(obs = 
                              (acyc_subset.longitude >= lon_min) & 
                              (acyc_subset.longitude <= lon_max) 
                             )<jupyter_output><empty_output><jupyter_text>Generate eddies geometries from the effective contour data in the META3.1exp database<jupyter_code># Cyclonic eddies
cyc_shape_lat = cyc_filter.effective_contour_latitude.data
cyc_shape_lon = (cyc_filter.effective_contour_longitude.data + 180) % 360 - 180
cyc_geometry = []

for x in range(len(cyc_shape_lat)):
    cyc_geometry.append(
        Polygon(
            zip(cyc_shape_lon[x], cyc_shape_lat[x])
        ))
# Anticyclonic eddies
acyc_shape_lat = acyc_filter.effective_contour_latitude.data
acyc_shape_lon = (acyc_filter.effective_contour_longitude.data + 180) % 360 - 180
acyc_geometry = []

for x in range(acyc_shape_lat.shape[0]):
    acyc_geometry.append(
        Polygon(
            zip(acyc_shape_lon[x], acyc_shape_lat[x])
        ))<jupyter_output><empty_output><jupyter_text>Create the geodataframes (gdf) for the eddies<jupyter_code># Cyclonic eddies
cyc_dic = {"track": cyc_filter.track,
       "time": cyc_filter.time,
       "longitude": cyc_filter.longitude,
       "latitude": cyc_filter.latitude,
       "area": cyc_filter.effective_area / 1000000, # 1000000 convert m² to km²
       "radius": cyc_filter.effective_radius / 1000, # 1000 convert m to km
       "geometry": cyc_geometry
      }

cyc_gdf = gpd.GeoDataFrame(cyc_dic, crs="EPSG:4326")
# Anticyclonic eddies
acyc_dic = {"track": acyc_filter.track,
        "time": acyc_filter.time,
        "longitude": acyc_filter.longitude,
        "latitude": acyc_filter.latitude,
        "area": acyc_filter.effective_area / 1000000,
        "radius": acyc_filter.effective_radius / 1000,
        "geometry": acyc_geometry
       }

acyc_gdf = gpd.GeoDataFrame(acyc_dic, crs = "EPSG:4326")<jupyter_output><empty_output><jupyter_text>Load the hydrographic stations position data for the TAI cruise and transform into a gdf<jupyter_code># Read the stations
positions = pd.read_csv("Cruise.csv")

# Transform the time to the same format used in the eddies
positions.Date = [np.datetime64(date, "D") for date in positions.Date]

# Create the geometries from the point coordinates
points_geometry = []
for x in range(positions.shape[0]):
    points_geometry.append(Point(positions.Longitude[x], positions.Latitude[x]))

# Create the geodataframe
stations_gdf = gpd.GeoDataFrame(positions, geometry = points_geometry, crs = "EPSG:4326")<jupyter_output><empty_output><jupyter_text>Make a spatial join between the eddies gdf and the cruise gdf and filter by date<jupyter_code># Cyclonic eddies
# Make a inner join, selecting only those station where is an intersection with an eddie
cyc_join = stations_gdf.sjoin(cyc_gdf, how = "inner")

# Filter by date, this way only the stations inside a eddie are selected
cyc_stations = cyc_join[cyc_join.Date == cyc_join.time]

# Generate a list of the sampled eddies tracks
cyc_unique = np.unique(cyc_stations.track)

# Extract from the META3.1 database the tracks information and create a list
cyc_tracks = [cyc.sel(obs = (cyc.track == x)) for x in cyc_unique]
# Anticyclonic eddies
# Make a inner join, selecting only those station where is an intersection with an eddie
acyc_join = stations_gdf.sjoin(acyc_gdf, how = "inner")

# Filter by date, this way only the stations inside a eddie are selected
acyc_stations = acyc_join[acyc_join.Date == acyc_join.time]

# Generate a list of the sampled eddies tracks
acyc_unique = np.unique(acyc_stations.track)

# Extract from the META3.1 database the tracks information and create a list
acyc_tracks = [acyc.sel(obs = (acyc.track == x)) for x in acyc_unique]<jupyter_output><empty_output><jupyter_text>*** **Part 2: Results** Stations from the TAI cruise inside an cyclonic eddie<jupyter_code>cyc_stations<jupyter_output><empty_output><jupyter_text>Stations from the TAI cruise inside an anticyclonic eddie<jupyter_code>acyc_stations<jupyter_output><empty_output><jupyter_text>Eddies sampled in the TAI cruise<jupyter_code># Create the figure, define the projection and the geographic extention
fig = plt.figure(figsize = (10,10), dpi = 300)
ax = plt.axes(projection = ccrs.PlateCarree())
ax.set_extent((-51, 21, -36, -14))

# Load the coastline and the grids in the map 
ax.coastlines(resolution = '50m')
ax.add_wms(
    wms="https://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?",
    layers=["GEBCO_LATEST"],
    alpha = 0.45
    )
gl = ax.gridlines(draw_labels = True, alpha = 0.1, color = "white")
gl.top_labels = False
gl.right_labels = False

# Add the sampled cyclonic eddies
cyc_filtered = cyc_gdf.loc[cyc_stations.index_right, :].drop_duplicates("track")
ax.add_geometries(cyc_filtered.geometry, crs = ccrs.PlateCarree(), facecolor='b', alpha=0.7)

# Add the sampled anticyclonic eddies
acyc_filtered = acyc_gdf.loc[acyc_stations.index_right, :].drop_duplicates("track")
ax.add_geometries(acyc_filtered.geometry, crs = ccrs.PlateCarree(), facecolor='r', alpha=0.7)

# Add the cruise stations
plt.plot(stations_gdf.Longitude, 
         stations_gdf.Latitude, 
         color = "black", 
         marker = '.', 
         transform = ccrs.PlateCarree(), 
         linewidth= 0,
         alpha = 0.8,
         markeredgewidth = 0)
 
plt.show()<jupyter_output><empty_output><jupyter_text>Complete tracks of the eddies sampled in the TAI cruise<jupyter_code># Create the figure, define the projection and the geographic extention
fig = plt.figure(figsize = (10,10), dpi = 300)
ax = plt.axes(projection = ccrs.PlateCarree())
ax.set_extent((-51, 21, -36, -14))

# Load the coastline and the grids in the map 
ax.coastlines(resolution = '50m')
ax.add_wms(
    wms="https://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?",
    layers=["GEBCO_LATEST"],
    alpha = 0.45
    )
gl = ax.gridlines(draw_labels = True, alpha = 0.15, color = "white")
gl.top_labels = False
gl.right_labels = False

# Create custom colors for the tracks
reds = ['#96222C', '#A62630', '#B72A35', '#C82D3A', '#D23744', '#D54854', '#D95963'] # OrRd Colorbrewer
blues = ["#20285B", "#252E6A", "#2A3579", "#303B88", "#354097", "#3A47A6", "#404DB5", "#404DB5", "#5966C5", "#6873CA"] # PuBu Colorbrewer

# Add the cyclonic tracks
for n in range(len(cyc_unique)):
    plt.plot(cyc_tracks[n].longitude, 
             cyc_tracks[n].latitude,
             color = list(reversed(blues))[n], 
             linewidth = 1.5, 
             transform = ccrs.PlateCarree()
            )
    plt.plot(cyc_tracks[n].longitude[0], 
             cyc_tracks[n].latitude[0],
             color = list(reversed(blues))[n], 
             marker = 'o', 
             transform = ccrs.PlateCarree()
            )
    plt.text(cyc_tracks[n].longitude[0] + 3.5, 
             cyc_tracks[n].latitude[0], 
             cyc_unique[n],
             fontsize = "xx-small",
             color = "b",
             horizontalalignment = 'right',
             transform = ccrs.PlateCarree()
            )

# Add the anticyclonic tracks
for n in range(len(acyc_unique)):
    plt.plot(acyc_tracks[n].longitude, 
             acyc_tracks[n].latitude, 
             color = list(reversed(reds))[n], 
             linewidth = 1.5, 
             transform = ccrs.PlateCarree()
            )
    plt.plot(acyc_tracks[n].longitude[0], 
             acyc_tracks[n].latitude[0], 
             color = list(reversed(reds))[n], 
             marker = 'o', 
             transform = ccrs.PlateCarree()
            )
    plt.text(acyc_tracks[n].longitude[0] + 3.5, 
             acyc_tracks[n].latitude[0], 
             acyc_unique[n],
             fontsize = "xx-small",
             color = "r",
             horizontalalignment = 'right',
             transform = ccrs.PlateCarree()
            )

# Show the figure
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11304178055447091, 'fraction_numerical': 0.030066380320187425, 'mean_word_length': 3.284818067754078, 'pattern_counts': {'":': 14, '<': 53, '<?xml version=': 0, '>': 55, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '574475', 'license': 'Apache-2.0', 'n_tokens_mistral': 3808, 'n_tokens_neox': 3458, 'n_words': 801, 'path': 'eddies_identification.ipynb', 'repo_name': 'andresawa/eddies-identification', 'size': '10263'}