Commit 746c5894 authored by David Vincent's avatar David Vincent Committed by David Vincent
Browse files

fix error in slimPre.rotateSphere

coordinates[0] are coordinates of point 0
coordinates[:,0] is an arry of longitudes of all points
creat slimPre.rotate3D which allow for rotating a vector from lon/lat coordinate system to the mesh coordinate system on a manifold

 fix function rotate in class Coordinate_system:
    the angle is now correctly compute and the transformation matrix is modified accordingly (the old function used an angle improperly computed)
parent fc9e0d46
......@@ -200,6 +200,7 @@ class Coordinate_system :
data projection system.
"""
self._data_proj = data_proj
self._region = region
if region._empty:
return
......@@ -226,14 +227,14 @@ class Coordinate_system :
mesh_proj2 = "+proj=latlong +ellps=WGS84"
_mesh_proj = dgpy.pj_init_plus(mesh_proj2)
_data_proj = dgpy.pj_init_plus(data_proj)
xyz = region.coordinates
if mesh_proj :
self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, xyz)
epsVect = np.zeros(xyz.shape)
epsVect[:,0] = np.hypot(xyz[1,0]-xyz[0,0],xyz[1,1]-xyz[0,1]) / 1000.
xyz0 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz-epsVect)
xyz1 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz+epsVect)
coord = self.coordinates.copy()
epsVect[:,0] = np.hypot(coord[1,0]-coord[0,0],coord[1,1]-coord[0,1]) / 1000.
xyz0 = dgpy.pjTransform(_data_proj, _mesh_proj, coord-epsVect)
xyz1 = dgpy.pjTransform(_data_proj, _mesh_proj, coord+epsVect)
if dgpy.pj_is_latlong(_data_proj) :
xyz1[:,0] = np.where( xyz1[:,0]<xyz0[:,0], xyz1[:,0]+2*np.pi, xyz1[:,0])
self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
......@@ -243,12 +244,37 @@ class Coordinate_system :
lonlat[:,1] = np.arcsin(xyz[:,2] / np.sqrt(xyz[:,0]*xyz[:,0]+xyz[:,1]*xyz[:,1]+xyz[:,2]*xyz[:,2]) )
lonlat[:,2] = np.sqrt(xyz[:,0]*xyz[:,0]+xyz[:,1]*xyz[:,1]+xyz[:,2]*xyz[:,2])
self.coordinates = dgpy.pjTransform(_mesh_proj, _data_proj, lonlat)
epsVect = np.zeros(xyz.shape)
epsVect[:,0] = 0.01
xyz0 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz-epsVect)
xyz1 = dgpy.pjTransform(_mesh_proj, _data_proj, xyz+epsVect)
self._angle = np.arctan2(xyz1[:,1]-xyz0[:,1], xyz1[:,0]-xyz0[:,0])
dlon = 1e-6*np.pi/180.
dlat = 1e-6*np.pi/180.
dr = 1e-6*np.pi/180.
epsVect = np.zeros(xyz.shape)
alpha = np.zeros(shape=(xyz.shape[0], xyz.shape[1],3))
epsVect[:,0] = dlon
epsVect[:,1] = dlat
epsVect[:,2] = dr
coord0 = self.coordinates[:,:]
coord1 = coord0.copy()
coord2 = coord0.copy()
coord3 = coord0.copy()
epsVect[coord1[:,0]+dlon>np.pi,0] = dlon-2.*np.pi
epsVect[coord1[:,1]+dlat>0.5*np.pi,1] = -dlat
coord1[:,0] = coord1[:,0] + epsVect[:,0]
coord2[:,1] = coord2[:,1] + epsVect[:,1]
coord3[:,2] = coord3[:,2] + epsVect[:,2]
xyz0 = dgpy.pjTransform(_data_proj, _mesh_proj, coord0)
xyz1 = dgpy.pjTransform(_data_proj, _mesh_proj, coord1)
xyz2 = dgpy.pjTransform(_data_proj, _mesh_proj, coord2)
xyz3 = dgpy.pjTransform(_data_proj, _mesh_proj, coord3)
for i in range(3):
alpha[:,i,0] = (xyz1[:,i]-xyz0[:,i])/epsVect[:,0]
alpha[:,i,1] = (xyz2[:,i]-xyz0[:,i])/epsVect[:,1]
alpha[:,i,2] = (xyz3[:,i]-xyz0[:,i])/epsVect[:,2]
alphaNorm = np.linalg.norm(alpha, axis = 1)
for i in range(3):
alpha[:,i,0] = alpha[:,i,0]/alphaNorm[:,0]
alpha[:,i,1] = alpha[:,i,1]/alphaNorm[:,1]
alpha[:,i,2] = alpha[:,i,2]/alphaNorm[:,2]
self._alpha = alpha
if _mesh_proj:
dgpy.pj_free(_mesh_proj)
if _data_proj:
......@@ -264,29 +290,47 @@ class Coordinate_system :
* v
second component of the vector (in the data coordinate system)
"""
u2 = u * np.cos(self._angle) + v * np.sin(self._angle)
v2 = -u * np.sin(self._angle) + v * np.cos(self._angle)
#u2 = u * np.cos(self._angle) - v * np.sin(self._angle)
#v2 = u * np.sin(self._angle) + v * np.cos(self._angle)
u2 = self._alpha[:,0,0]*u + self._alpha[:,0,1]*v
v2 = self._alpha[:,1,0]*u + self._alpha[:,1,1]*v
return (u2, v2)
def rotate3D(self, F_lon, F_colat, F_r):
"""Return ...
def rotateSphere(self, F_lon, F_colat, F_r):
"""Return the vector in the carthesian field of coordinates
keyword arguments:
* F_lon
this var ...
component along the longitude direction
* F_colat
this var ...
component along the colatitude direction
* F_r
this var ...
component along the radius direction
"""
lon = self.coordinates[0]
colat = 0.5*np.pi - self.coordinates[1]
lon = self.coordinates[:,0]
colat = 0.5*np.pi - self.coordinates[:,1]
Fx = np.sin(colat)*np.cos(lon)*F_r + np.cos(colat)*np.cos(lon)*F_colat - np.sin(lon)*F_lon
Fy = np.sin(colat)*np.sin(lon)*F_r + np.cos(colat)*np.sin(lon)*F_colat + np.cos(colat)*F_lon
Fy = np.sin(colat)*np.sin(lon)*F_r + np.cos(colat)*np.sin(lon)*F_colat + np.cos(lon)*F_lon
Fz = np.cos(colat)*F_r - np.sin(colat)*F_colat
return (Fx,Fy,Fz)
def rotate3D(self, u, v, w):
"""Return the rotated vector data array (u,v). Returned vector will be oriented in the carthesian system
keyword arguments:
* u
first component of the vector (in the data coordinate system)
* v
second component of the vector (in the data coordinate system)
* w
third component of the vector (in the data coordinate system)
"""
u2 = self._alpha[:,0,0]*u + self._alpha[:,0,1]*v + self._alpha[:,0,2]*w
v2 = self._alpha[:,1,0]*u + self._alpha[:,1,1]*v + self._alpha[:,1,2]*w
w2 = self._alpha[:,2,0]*u + self._alpha[:,2,1]*v + self._alpha[:,2,2]*w
return (u2, v2, w2)
class Time :
"""Create the slimPre Time"""
......
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