#include #include #include #include #include "proj.h" #include "grb2.h" #include "wgrib2.h" #include "fnlist.h" /* gctpc_aux.c interface routines to the gctpc library 2/2012 Public Domain Wesley Ebisuzaki gctpc_get_latlon: fill grid with lat/lon values mercator polar stereographic lambert conformal Albers equal area 10/2015 lambert azimuthal equal area */ /* M_PI, M_PI_2, M_PI_4, and M_SQRT2 are not ANSI C but are commonly defined */ /* values from GNU C library version of math.h copyright Free Software Foundation, Inc. */ #ifndef M_PI #define M_PI 3.14159265358979323846 /* pi */ #endif #ifndef M_PI_2 #define M_PI_2 1.57079632679489661923 /* pi/2 */ #endif #ifndef M_PI_4 #define M_PI_4 0.78539816339744830962 /* pi/4 */ #endif #ifndef M_SQRT2 #define M_SQRT2 1.41421356237309504880 /* sqrt(2) */ #endif /* * given vectors of lon-lat * return i-j vectors * * step 1: find (x,y) of 1st point * step 2: find (x,y) of input point * * input: * **sec = grib sections of grid * **grid_lon = longitues of grid * **grid_lat = latitudes of grid * n = number of point to convert * *lon = longitudes of points to convert to (i,j) * *lat = longitudes of points to convert to (i,j) * * output: * x[n] coordinates on grid * y[n] * * assumption: grid_lon, grid_lat is in wesn order * * to use: * * setup * int gctpc_ll2xy_init(unsigned char **sec, double *grid_lon, double *grid_lat); * * get (x,y) * int gctpc_ll2xy(int n, double *lon, double *lat, double *x, double *y); * get nearest neighbor * int gctpc_ll2i(int n, double *lon, double *lat, unsigned int *i); * */ static long int (*forward_fn)(); static double dx, dy, inv_dx, inv_dy, x_0, y_0, x00, xN; static unsigned int gdt; static int nx, ny; int gctpc_ll2xy_init(unsigned char **sec, double *grid_lon, double *grid_lat) { unsigned char *gds; double r_maj; /* major axis */ double r_min; /* minor axis */ double lat1; /* first standard parallel */ double lat2; /* second standard parallel */ double c_lon; /* center longitude */ double c_lat; /* center latitude */ double false_east; /* x offset in meters */ double false_north; int nres, nscan; unsigned int npnts; long long_i; double rlon, rlat; gdt = code_table_3_1(sec); gds = sec[3]; /* only process certain grids */ forward_fn = NULL; if (grid_lon == NULL || grid_lat == NULL) return 1; if (gdt != 0 && !(gdt == 10 && GDS_Mercator_ori_angle(gds) == 0.0) && gdt != 20 && gdt != 30 && gdt != 140) return 1; get_nxny(sec, &nx, &ny, &npnts, &nres, &nscan); if (nx == -1 || ny == -1 || nx*ny != npnts) return 1; /* only process certain grids */ x_0 = y_0 = x00 = xN = inv_dx = inv_dy = 0.0; if (gdt == 0) { /* lat-lon grid */ dx = grid_lon[1] - grid_lon[0]; dy = grid_lat[nx] - grid_lat[0]; inv_dx = 1.0 / dx; inv_dy = 1.0 / dy; x_0 = grid_lon[0]; x00 = grid_lon[0] - 0.5*dx; xN = grid_lon[nx-1] + 0.5*dx; y_0 = grid_lat[0]; return 0; } else if (gdt == 10 && (GDS_Mercator_ori_angle(gds) == 0.0) ) { // mercator no rotation /* get earth axis */ axes_earth(sec, &r_maj, &r_min); dx = GDS_Mercator_dx(gds); dy = GDS_Mercator_dy(gds); inv_dx = 1.0 / dx; inv_dy = 1.0 / dy; /* central point */ c_lon = 0.0; c_lat = GDS_Mercator_latD(gds) * (M_PI/180.0); /* find the eastling and northing of of the 1st grid point */ false_east = false_north = 0.0; long_i = merforint(r_maj,r_min,c_lon,c_lat,false_east,false_north); if (long_i) fatal_error_i("gctpc_ll2xy_init merforint: return %ld", long_i); rlat = GDS_Mercator_lat1(gds) * (M_PI/180.0); rlon = GDS_Mercator_lon1(gds) * (M_PI/180.0); long_i = merfor(rlon, rlat, &x_0, &y_0); x00 = x_0 - 0.5*dx; xN = x_0 + (nx-0.5)*dx; forward_fn = &merfor; } else if (gdt == 20) { // polar stereographic /* get earth axis */ axes_earth(sec, &r_maj, &r_min); dy = GDS_Polar_dy(gds); dx = GDS_Polar_dx(gds); inv_dx = 1.0 / dx; inv_dy = 1.0 / dy; /* central point */ c_lon = GDS_Polar_lov(gds) * (M_PI/180.0); c_lat = GDS_Polar_lad(gds) * (M_PI/180.0); /* find the eastling and northing of of the 1st grid point */ false_east = false_north = 0.0; long_i = psforint(r_maj,r_min,c_lon,c_lat,false_east,false_north); if (long_i) fatal_error_i("gctpc_ll2xy_init psforint: return %ld", long_i); rlon = grid_lon[0] * (M_PI/180.0); rlat = grid_lat[0] * (M_PI/180.0); long_i = psfor(rlon, rlat, &x_0, &y_0); x00 = x_0 - 0.5*dx; xN = x_0 + (nx-0.5)*dx; forward_fn = &psfor; } else if (gdt == 30) { // lambert conformal conic /* get earth axis */ axes_earth(sec, &r_maj, &r_min); dx = GDS_Lambert_dx(gds); dy = GDS_Lambert_dy(gds); inv_dx = 1.0 / dx; inv_dy = 1.0 / dy; /* latitudes of tangent/intersection */ lat1 = GDS_Lambert_Latin1(gds) * (M_PI/180.0); lat2 = GDS_Lambert_Latin2(gds) * (M_PI/180.0); /* central point */ c_lon = GDS_Lambert_Lov(gds) * (M_PI/180.0); c_lat = GDS_Lambert_LatD(gds) * (M_PI/180.0); /* find the eastling and northing of of the 1st grid point */ false_east = false_north = 0.0; long_i = lamccforint(r_maj,r_min,lat1,lat2,c_lon,c_lat,false_east,false_north); if (long_i) fatal_error_i("gctpc_ll2xy_init lamccforint: return %ld", long_i); rlon = grid_lon[0] * (M_PI/180.0); rlat = grid_lat[0] * (M_PI/180.0); long_i = lamccfor(rlon, rlat, &x_0, &y_0); x00 = x_0 - 0.5*dx; xN = x_0 + (nx-0.5)*dx; forward_fn = &lamccfor; } else if (gdt == 140) { // lambert azimuthal equal area /* get earth axis */ axes_earth(sec, &r_maj, &r_min); r_maj = 0.5 * (r_maj + r_min); dx = GDS_Lambert_Az_dx(gds); dy = GDS_Lambert_Az_dy(gds); inv_dx = 1.0 / dx; inv_dy = 1.0 / dy; /* central point */ c_lon = GDS_Lambert_Az_Cen_Lon(gds) * (M_PI/180.0); c_lat = GDS_Lambert_Az_Std_Par(gds) * (M_PI/180.0); false_east = false_north = 0.0; long_i = lamazforint(r_maj,c_lon,c_lat,false_east,false_north); if (long_i) fatal_error_i("gctpc_ll2xy_init lamazforint: return %ld", long_i); rlon = grid_lon[0] * (M_PI/180.0); rlat = grid_lat[0] * (M_PI/180.0); long_i = lamazfor(rlon, rlat, &x_0, &y_0); x00 = x_0 - 0.5*dx; xN = x_0 + (nx-0.5)*dx; forward_fn = &lamazfor; } return forward_fn != NULL ? 0 : 1; } int gctpc_ll2xy(int n, double *lon, double *lat, double *x, double *y) { int i; double rlon, rlat; if (gdt == 0) { // lat-lon // #pragma omp parallel for schedule(static) private(i,rlon,rlat) for (i = 0; i < n; i++) { rlon = lon[i]; if (rlon > xN) rlon -= 360.0; if (rlon < x00) rlon += 360.0; rlat = lat[i]; x[i] = (rlon - x_0) * inv_dx; y[i] = (rlat - y_0) * inv_dy; } return 0; } if (forward_fn == NULL) return 1; // #pragma omp parallel for private(i,rlon,rlat) for (i = 0; i < n; i++) { rlon = lon[i]; if (rlon > xN) rlon -= 360.0; if (rlon < x00) rlon += 360.0; rlon *= (M_PI/180.0); rlat = lat[i] * (M_PI/180.0); forward_fn(rlon, rlat, x+i, y+i); x[i] = (x[i] - x_0)*inv_dx; y[i] = (y[i] - y_0)*inv_dy; } return 0; } /* iptr[] == 0 for out of bounds */ int gctpc_ll2i(int n, double *lon, double *lat, unsigned int *ipnt) { int i; unsigned int ix, iy; double rlon, rlat, x, y; if (gdt == 0) { // lat-lon // #pragma omp parallel for schedule(static) private(i,rlon,rlat,ix,iy, x, y) for (i = 0; i < n; i++) { rlon = lon[i]; if (rlon > xN) rlon -= 360.0; if (rlon < x00) rlon += 360.0; rlat = lat[i]; ix = x = floor((rlon - x_0) * inv_dx + 0.5); iy = y = floor((rlat - y_0) * inv_dy + 0.5); if (x < 0 || x >= nx || y < 0 || y >= ny) { ipnt[i] = 0; } else { ipnt[i] = ix + nx*iy + 1; } } return 0; } if (forward_fn == NULL) return 1; // #pragma omp parallel for schedule(static) private(i,rlon,rlat,ix,iy, x, y) for (i = 0; i < n; i++) { rlon = lon[i]; if (rlon > xN) rlon -= 360.0; if (rlon < x00) rlon += 360.0; rlon *= (M_PI/180.0); rlat = lat[i] * (M_PI/180.0); forward_fn(rlon, rlat, &x, &y); ix = x = floor((x - x_0)*inv_dx + 0.5); iy = y = floor((y - y_0)*inv_dy + 0.5); if (x < 0 || x >= nx || y < 0 || y >= ny) { ipnt[i] = 0; } else { ipnt[i] = ix + nx*iy + 1; } } return 0; }