Locating shearlines or deformation zones in a map

Question

I am looking for a way to detect shearlines or confluence regions such as that shown in the figure below (red dashed line).

I was able to derive the stretching, shear, and total deformations.

The shearline coincides with the region of maximum stretching deformation as in the figure below (red dashed line):

Currently, I am drawing the shearline manually. But I am wondering if there is a way to automatically detect this line or the deformation zone.

I am thinking of two ways

1. Looking at the dynamic variables (i.e., stretching deformation)

2. Using Python's pattern recognition function to detect this kind of pattern.

Problems

Any suggestions on how I can detect this automatically. I don't know how to implement the above methods. I am not sure what criteria to set in order to draw the shearline.

I am also new to Python so I don't know how to implement the second one.

I have wind data and calculated deformations in netcdf format. I uploaded the files to the following links:

UWND

VWND

Link to calculated deformation parameters

The deformation parameters are calculated this way:

where A, B, and V are the shearing, stretching, and total deformation parameters.

In the deform.nc file, shear_sph(shearing deformation; A in the equation above); stretch_sph(stretching deformation; B in the equation B); deform_sph(total deformation; V in the equation above).

I'll appreciate any help about this matter.

SUMMARY

I am looking for the NCL Shear, stretch, deform functionality in Python. Are there any existing software packages or do I have to write my own ?

Answer 1

I am not able to give full actual code in Python but I hope I can point you in the right direction if I were to go about calculating diffluence and confluence.

One can start with making a streamlines plot using matplotlib first. Once you have that you need to get the actual streamlines data as described in this link - https://matplotlib.org/api/_as_gen/matplotlib.pyplot.streamplot.html and that data structure is given as LineCollection.

For each streamline in that LineCollection you will need to calculate tangential and normal vectors based on the idea of a localized Cartesian system with a origin centered on somewhere along the streamline.

Once you have the tangential and normal vectors you will need to use this equation i.e. the divergence in natural coordinates

$$ \nabla . v_h = (e_n \frac {\partial e_s}{\partial n} * v_h) + \frac {\partial v_h}{\partial s}$$

Using a combination of Frenet Serret formula this can be shown to be equal to

$$\frac {|v_h|}{R_n} + \frac {\partial v_h}{\partial s} $$

From your perspective you only want the first term on that equation to calculate diffluence and confluence. So $R_n$ is the radius of curvature of the normal vector (which you already calculated before). So in order to calculate the radius of curvature you can use this python code for least squares fitting (using the approach of least squares - obtain the ellipse and then calculate the radius). If the radius is positive i..e $R_n >0 $ then you have an area of diffluence and negative then confluence.

UPDATE

If you are looking for the exact NCL shear stretch deform in Python what you can do is to use this package in python pyspharm. The data has to be global and then you need to orient the grid in a south to north direction. Once you have that you can use that package to calculate the shear, stretch and deformation parameters by using the pyspharm's getgrad() method to calculate the vector gradient of the wind vectors. If you pass in the components of the wind vectors(both u and v) to getgrad() you will get the zonal and meridional components of the gradient of the wind vectors. Adding and subtracting them will get you the shear and stretch.

Sample code

 from spharm import Spharmt

 uspectral = Spharmt.grdtospec(u,ntruncation)
 vspectral = Spharmt.grdtospec(v,ntruncation)
 dudx,dudy = Spharmt.getgrad(uspectral)
 dvdx,dvdy = Spharmt.getgrad(vspectral)
 #shear and stretch
 shear = dvdx + dudy
 stretch = dudx - dvdy