Commit f61d2fbf authored by CEVAER's avatar CEVAER
Browse files

Fixed wrong lon/lat grids

parent 1e3a0431
......@@ -7,6 +7,7 @@ import os
import logging
logger = logging.getLogger(__name__)
def research_contour(mask, lon, lat):
idx0, idx1 = contour_indices(mask, lon)
lon_lwind = lon[idx0, idx1]
......@@ -37,7 +38,7 @@ def save_netcdf(sat_file, path, lon_center, lat_center, lons_eye_shape, lats_eye
ds_polar_north.attrs["High wind area barycenter"] = Point(lon_eye_hwind, lat_eye_hwind).wkt
ds_polar_north.attrs["Interpolated ATCF track point"] = Point(lon_eye_fg, lat_eye_fg).wkt
ds_polar_north.attrs["Storm name"] = storm_name
ds_polar_north.attrs["ATCF VMAX"] = float(atcf_vmax)
ds_polar_north.attrs["Track VMAX"] = float(atcf_vmax) # m/s
lon_lwind, lat_lwind = research_contour(msk_lwind, lon, lat)
lwind_contour_coords = []
......@@ -53,6 +54,12 @@ def save_netcdf(sat_file, path, lon_center, lat_center, lons_eye_shape, lats_eye
hwind_contour = LineString(hwind_contour_coords)
ds_polar_north.attrs["High wind research area"] = hwind_contour.wkt
outfile = os.path.join(path, f"{os.path.basename(sat_file).replace('.nc', '').replace('cm', 'ca')}.nc")
# Storm id
# Source track file
# source sat file
# storm id nom du fichier
# vmax m/s
outfile = os.path.join(path, f"{os.path.basename(sat_file_sw).replace('.nc', '').replace('cm', 'ca')}.nc")
ds_polar_north.to_netcdf(path=outfile)
logger.info(f"Wrote netCDF product to {outfile}")
......@@ -2,7 +2,7 @@ import cv2
import numpy as np
import xarray as xr
from pyproj import Geod
import time
def find_nearest_2d(lon, lat, lonc, latc):
abslat = np.abs(lat - latc)
......@@ -14,8 +14,29 @@ def find_nearest_2d(lon, lat, lonc, latc):
return xloc, yloc
def pixel_distance():
pass
def pixel_distance(ds, xloc, yloc):
lon1 = ds.lon.values[xloc, yloc]
lat1 = ds.lat.values[xloc, yloc]
lon2 = ds.lon.values[xloc + 1, yloc]
lat2 = ds.lat.values[xloc + 1, yloc]
geod = Geod(ellps='WGS84')
faz, baz, dist = geod.inv(lon1, lat1, lon2, lat2)
return dist
def custom_meshgrid(theta, rad, data_shape):
st = time.time()
coord_theta = np.zeros(data_shape)
coord_rad = np.zeros(data_shape)
for i in range(0, data_shape[0]):
for j in range(0, data_shape[1]):
coord_theta[i, j] = theta[i]
coord_rad[i, j] = rad[j]
print("Time taken ", time.time() - st)
return coord_theta, coord_rad
def to_polar(ds, lon_center, lat_center):
......@@ -25,18 +46,17 @@ def to_polar(ds, lon_center, lat_center):
# Find indexes in data array of the lon/lat center
xloc, yloc = find_nearest_2d(ds.lon, ds.lat, lon_center, lat_center)
pix_dist = pixel_distance(ds, xloc, yloc)
print(pix_dist)
# Coordinates array
theta = np.linspace(0, 360, num=data_shape[0])
rads = np.linspace(0, max_radius, num=data_shape[1])
rads_meter = rads * pix_dist
print(rads_meter.min(), rads_meter.max())
coord_theta = np.zeros(data_shape)
coord_rad = np.zeros(data_shape)
for i in range(0, data_shape[0]):
for j in range(0, data_shape[1]):
coord_theta[i, j] = theta[i]
coord_rad[i, j] = rads[j]
#coord_theta, coord_rad = np.meshgrid(theta, rads_meter)
coord_theta, coord_rad = np.meshgrid(theta, rads)
coord_theta, coord_rad = custom_meshgrid(theta, rads_meter, data_shape)
geod = Geod(ellps='WGS84')
lons_polar, lats_polar, baz = geod.fwd(ds.lon.values * 0. + lon_center, ds.lat.values * 0. + lat_center,
......@@ -53,7 +73,7 @@ def to_polar(ds, lon_center, lat_center):
polar_ds = xr.Dataset(data_vars=vars_polar, coords={
"theta": (["theta"], theta),
"rad": (["rad"], rads)
"rad": (["rad"], rads_meter)
})
return polar_ds
......@@ -80,7 +100,7 @@ def to_north_ref(polar_ds):
polar_ds.theta.values[:] = np.roll(polar_ds.theta.values, rolling)
# Roll the whole ds back so that theta[0] = 0
polar_ds = polar_ds.roll(theta=-rolling)
polar_ds = polar_ds.roll(theta=-rolling, roll_coords=True)
return polar_ds
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment