Closed gallego2286 closed 1 year ago
To create shapefile you could use python module "shapefile", there is no routine in py eddy tracker to create shapefile. You will need to iterate on effective_contour_longitude /effective_contour_latitude or/and speed_contour_longitude/speed_contour_latitude in order to create shape for each eddy.
To create shapefile you could use python module "shapefile", there is no routine in py eddy tracker to create shapefile. You will need to iterate on effective_contour_longitude /effective_contour_latitude or/and speed_contour_longitude/speed_contour_latitude in order to create shape for each eddy.
Thank you very much for your timely and quick response. Best regards.
If you wrote a method to do this you are welcome to share it here
Of course, I hope this can help someone else if they find themselves in the same situation as I did a couple of days ago. Thank you very much for your help again.
# %%
import os
from datetime import datetime
from matplotlib import pyplot as plt
from py_eddy_tracker.dataset.grid import RegularGridDataset
from py_eddy_tracker.observations.observation import EddiesObservations
from netCDF4 import Dataset
import geopandas as gpd
from shapely import box, LineString
# %%
#path of the example images
directory = r'c:\Users\user\Desktop\py-eddy-tracker-master\src\py_eddy_tracker\data'
#sample file name
file_example = f'nrt_global_allsat_phy_l4_20190223_20190226.nc'
complete_path = os.path.join(directory, file_example)
# %%
from datetime import datetime
h = RegularGridDataset(complete_path, "longitude", "latitude")
h.bessel_high_filter("adt", 700, order=1)
date = datetime(2019, 2, 23)
a, c = h.eddy_identification(
"adt",
"ugos",
"vgos", # Variables used for identification
date, # Date of identification
0.002, # step between two isolines of detection (m)
pixel_limit=(5, 2000), # Min and max pixel count for valid contour
shape_error=55, # Error max (%) between ratio of circle fit and contour
)
# %%
fig = plt.figure(figsize=(15, 7))
ax = fig.add_axes([0.03, 0.03, 0.94, 0.94])
ax.set_title("Eddies detected -- Cyclonic(red) and Anticyclonic(blue)")
ax.set_ylim(-75, 75)
ax.set_xlim(0, 360)
ax.set_aspect("equal")
a.display(ax, color="b", linewidth=0.5)
c.display(ax, color="r", linewidth=0.5)
ax.grid()
#fig.savefig("share/png/eddies.png")
# %%
#saving path for processed images
stored_directory = r'c:\Users\user\Desktop\process'
#name of the example file
anticyclo = f'anticyclo_20190223_20190226.nc'
cyclonic = f'cyclonic_20190223_20190226.nc'
anticyclo_path = os.path.join(stored_directory, anticyclo)
cyclonic_path = os.path.join(stored_directory, cyclonic)
# %%
#Saving of processed files
with Dataset(anticyclo_path, "w") as h: a.to_netcdf(h)
with Dataset(cyclonic_path, "w") as h: c.to_netcdf(h)
# %%
# Extraction of the coordinates of the vertices
# of the detected cyclonic eddies
# Load cyclonic and anticyclonic files
Anticyclonic = EddiesObservations.load_file(anticyclo_path)
Cyclonic = EddiesObservations.load_file(cyclonic_path)
#Here you can see that I only do it with anticyclones.
lat_anti = Anticyclonic['contour_lat_e']
lon_anti = Anticyclonic['contour_lon_e']
# %%
# Create a list of polylines from coordinates
poli_anti = [LineString(zip(lon, lat))
for lat, lon in zip(lat_anti, lon_anti)]
# GeoDataFrame creation
gdf = gpd.GeoDataFrame({'geometry': poli_anti})
#Field creation to distinguish the anticyclones
gdf['eddies'] = 'anticyclo'
gdf.plot()
Dear Authors.
First I want to thank you for this very useful library that you have created for the community. Thank you very much for your efforts.
The question I asked, since I don't see that someone has asked it before in this section, is aimed at the result transformation obtained from the identification. I am a beginner and a question has arisen with the handling of this example: https://py-eddy-tracker.readthedocs.io/en/latest/grid_identification.html in the end I have been able to save the generated files here:
from netCDF import Dataset
with Dataset(date.strftime("share/Anticyclonic_%Y%m%d.nc"), "w") as h: a.tonetcdf(h). with Dataset(date.strftime("share/Cyclonic%Y%m%d.nc"), "w") as h: c.to_netcdf(h)
But from there, I have not been able to transform those files from .nc format to for example vector (shp) format, could you please give me a hand how can I do it. Thank you very much for your help.