core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Refactoring

This commit is contained in:
Shahram Najm 2020-06-14 20:10:44 +01:00
parent 49ca0c0167
commit 622fb168b0
6 changed files with 78 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/grib_api_prototypes.h
View File

@ -1048,6 +1048,8 @@ int grib_get_gaussian_latitudes(long trunc, double* lats);
int is_gaussian_global(double lat1, double lat2, double lon1, double lon2, long num_points_equator, const double* latitudes, double angular_precision);
void rotate(const double inlat, const double inlon, const double angleOfRot, const double southPoleLat, const double southPoleLon, double* outlat, double* outlon);
void unrotate(const double inlat, const double inlon, const double angleOfRot, const double southPoleLat, const double southPoleLon, double* outlat, double* outlon);
double geographic_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2);
double geographic_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2);
/* grib_handle.c */
grib_section* grib_section_create(grib_handle* h, grib_accessor* owner);
@ -1376,8 +1378,6 @@ int grib_nearest_find(grib_nearest* nearest, const grib_handle* h, double inlat,
int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args);
int grib_nearest_delete(grib_nearest* i);
void grib_binary_search(double xx[], const unsigned long n, double x, int* ju, int* jl);
double grib_nearest_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2);
double grib_nearest_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2);
int grib_nearest_find_multiple(const grib_handle* h, int is_lsm, const double* inlats, const double* inlons, long npoints, double* outlats, double* outlons, double* values, double* distances, int* indexes);
int grib_nearest_find_generic(grib_nearest* nearest, grib_handle* h, double inlat, double inlon, unsigned long flags,
const char* values_keyname, const char* radius_keyname,

70
src/grib_geography.c
View File

@ -4140,3 +4140,73 @@ void unrotate(const double inlat, const double inlon,
*outlat = ret_lat;
*outlon = ret_lon;
}
#define RADIAN(x) ((x)*acos(0.0) / 90.0)
/* radius is in km, angles in degrees */
/* Spherical Law of Cosines */
double geographic_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2)
{
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = lon1;
double rlon2 = lon2;
double a;
if (lat1 == lat2 && lon1 == lon2) {
return 0.0; /* the two points are identical */
}
if (rlon1 >= 360) rlon1 -= 360.0;
rlon1 = RADIAN(rlon1);
if (rlon2 >= 360) rlon2 -= 360.0;
rlon2 = RADIAN(rlon2);
a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
DebugAssert(a >= -1 && a <= 1);
return radius * acos(a);
}
/* major and minor axes in km, angles in degrees */
double geographic_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2)
{
/* Lambert's formula */
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = RADIAN(lon1);
double rlon2 = RADIAN(lon2);
double deltaLat = rlat2 - rlat1;
double deltaLon = rlon2 - rlon1;
double sinDlat = sin(deltaLat/2.0);
double sinDlon = sin(deltaLon/2.0);
double sin2Dlat = sinDlat*sinDlat;
double sin2Dlon = sinDlon*sinDlon;
double a = sin2Dlat + cos(rlat1) * cos(rlat2) * sin2Dlon;
double c = 2 * atan2(sqrt(a), sqrt(1.0-a));
double f = (major - minor)/major; /*flattening*/
double latr1 = atan( (1.0-f)*tan(rlat1) ); /*Reduced latitude1*/
double latr2 = atan( (1.0-f)*tan(rlat2) ); /*Reduced latitude2*/
double P = (latr1+latr2)/2;
double Q = (latr2-latr1)/2;
double sinP = sin(P);
double sin2P = sinP*sinP;
double cosQ = cos(Q);
double cos2Q = cosQ*cosQ;
double cosc2 = cos(c/2);
double cos2c2 = cosc2*cosc2;
double sinQ = sin(Q);
double sin2Q = sinQ*sinQ;
double cosP = cos(P);
double cos2P = cosP*cosP;
double sinc2 = sin(c/2);
double sin2c2 = sinc2*sinc2;
double X = (c - sin(c))* sin2P * cos2Q / cos2c2;
double Y = (c + sin(c))*sin2Q*cos2P/sin2c2;
double dist = major * (c - f*(X+Y)/2);
return dist;
}

72
src/grib_nearest.c
View File

@ -114,76 +114,6 @@ void grib_binary_search(double xx[], const unsigned long n, double x,
}
}
#define RADIAN(x) ((x)*acos(0.0) / 90.0)
/* radius is in km, angles in degrees */
/* Spherical Law of Cosines */
double grib_nearest_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2)
{
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = lon1;
double rlon2 = lon2;
double a;
if (lat1 == lat2 && lon1 == lon2) {
return 0.0; /* the two points are identical */
}
if (rlon1 >= 360) rlon1 -= 360.0;
rlon1 = RADIAN(rlon1);
if (rlon2 >= 360) rlon2 -= 360.0;
rlon2 = RADIAN(rlon2);
a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
DebugAssert(a >= -1 && a <= 1);
return radius * acos(a);
}
/* major and minor axes in km, angles in degrees */
double grib_nearest_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2)
{
/* Lambert's formula */
double rlat1 = RADIAN(lat1);
double rlat2 = RADIAN(lat2);
double rlon1 = RADIAN(lon1);
double rlon2 = RADIAN(lon2);
double deltaLat = rlat2 - rlat1;
double deltaLon = rlon2 - rlon1;
double sinDlat = sin(deltaLat/2.0);
double sinDlon = sin(deltaLon/2.0);
double sin2Dlat = sinDlat*sinDlat;
double sin2Dlon = sinDlon*sinDlon;
double a = sin2Dlat + cos(rlat1) * cos(rlat2) * sin2Dlon;
double c = 2 * atan2(sqrt(a), sqrt(1.0-a));
double f = (major - minor)/major; /*flattening*/
double latr1 = atan( (1.0-f)*tan(rlat1) ); /*Reduced latitude1*/
double latr2 = atan( (1.0-f)*tan(rlat2) ); /*Reduced latitude2*/
double P = (latr1+latr2)/2;
double Q = (latr2-latr1)/2;
double sinP = sin(P);
double sin2P = sinP*sinP;
double cosQ = cos(Q);
double cos2Q = cosQ*cosQ;
double cosc2 = cos(c/2);
double cos2c2 = cosc2*cosc2;
double sinQ = sin(Q);
double sin2Q = sinQ*sinQ;
double cosP = cos(P);
double cos2P = cosP*cosP;
double sinc2 = sin(c/2);
double sin2c2 = sinc2*sinc2;
double X = (c - sin(c))* sin2P * cos2Q / cos2c2;
double Y = (c + sin(c))*sin2Q*cos2P/sin2c2;
double dist = major * (c - f*(X+Y)/2);
return dist;
}
int grib_nearest_find_multiple(
const grib_handle* h, int is_lsm,
const double* inlats, const double* inlons, long npoints,
@ -447,7 +377,7 @@ int grib_nearest_find_generic(
/* Ignore latitudes too far from our point */
}
else {
double dist = grib_nearest_distance_spherical(radius, inlon, inlat, lon, lat);
double dist = geographic_distance_spherical(radius, inlon, inlat, lon, lat);
if (dist < min_dist)
min_dist = dist;
/*printf("Candidate: lat=%.5f lon=%.5f dist=%f Idx=%ld Val=%f\n",lat,lon,dist,the_index,the_value);*/

2
src/grib_nearest_class_latlon_reduced.c
View File

@ -346,7 +346,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
kk = 0;
for (jj = 0; jj < 2; jj++) {
for (ii = 0; ii < 2; ii++) {
self->distances[kk] = grib_nearest_distance_spherical(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->k[kk]], self->lats[self->j[jj]]);
kk++;
}

4
src/grib_nearest_class_reduced.c
View File

@ -409,7 +409,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
kk = 0;
for (jj = 0; jj < 2; jj++) {
for (ii = 0; ii < 2; ii++) {
self->distances[kk] = grib_nearest_distance_spherical(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->k[kk]], self->lats[self->j[jj]]);
kk++;
}
@ -548,7 +548,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
kk = 0;
for (ii = 0; ii < 2; ii++) {
for (jj = 0; jj < 2; jj++) {
self->distances[kk] = grib_nearest_distance_spherical(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->k[kk]], self->lats[self->j[jj]]);
kk++;
}

4
src/grib_nearest_class_regular.c
View File

@ -193,7 +193,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
for (ii=0;ii<2;ii++) {
for (jj=0;jj<2;jj++) {
self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1;
self->distances[kk]=grib_nearest_distance_spherical(radius,inlon,inlat,
self->distances[kk]=geographic_distance_spherical(radius,inlon,inlat,
self->lons[self->i[ii]],self->lats[self->j[jj]]);
kk++;
}
@ -409,7 +409,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
for (jj = 0; jj < 2; jj++) {
for (ii = 0; ii < 2; ii++) {
self->k[kk] = self->i[ii] + self->lons_count * self->j[jj];
self->distances[kk] = grib_nearest_distance_spherical(radius, inlon, inlat,
self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
self->lons[self->i[ii]], self->lats[self->j[jj]]);
kk++;
}