#include #include #include #include "mosaic_util.h" #include "create_xgrid.h" #include "constant.h" #define AREA_RATIO_THRESH (1.e-6) #define MASK_THRESH (0.5) #define EPSLN8 (1.e-8) #define EPSLN30 (1.0e-30) #define EPSLN10 (1.0e-10) #define R2D (180/M_PI) double grid_box_radius(const double *x, const double *y, const double *z, int n); double dist_between_boxes(const double *x1, const double *y1, const double *z1, int n1, const double *x2, const double *y2, const double *z2, int n2); int inside_edge(double x0, double y0, double x1, double y1, double x, double y); int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double *q3, double *intersect, double *u_a, double *u_q, int *inbound); /******************************************************************************* int get_maxxgrid return constants MAXXGRID. *******************************************************************************/ int get_maxxgrid(void) { return MAXXGRID; } int get_maxxgrid_(void) { return get_maxxgrid(); } /******************************************************************************* void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, const double *area) return the grid area. *******************************************************************************/ #ifndef __AIX void get_grid_area_(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area) { get_grid_area(nlon, nlat, lon, lat, area); } #endif void get_grid_area(const int *nlon, const int *nlat, const double *lon, const double *lat, double *area) { int nx, ny, nxp, i, j, n_in; double x_in[20], y_in[20]; nx = *nlon; ny = *nlat; nxp = nx + 1; for(j=0; j 1) get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in); else get_grid_area_no_adjust(nlon_in, nlat_in, tmpx, tmpy, area_in); get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out); free(tmpx); free(tmpy); for(j1=0; j1 MASK_THRESH ) { ll_lon = lon_in[i1]; ll_lat = lat_in[j1]; ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1]; for(j2=0; j2=ur_lat) && (y_in[1]>=ur_lat) && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_out[j2*nx2p+i2]; x_in[1] = lon_out[j2*nx2p+i2+1]; x_in[2] = lon_out[(j2+1)*nx2p+i2+1]; x_in[3] = lon_out[(j2+1)*nx2p+i2]; n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { Xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if( Xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = Xarea; i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area_in); free(area_out); return nxgrid; }; /* create_xgrid_1dx2d_order1 */ /******************************************************************************** void create_xgrid_1dx2d_order2 This routine generate exchange grids between two grids for the second order conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_in,lat_in are 1-D grid bounds, lon_out,lat_out are geographic grid location of grid cell bounds. ********************************************************************************/ int create_xgrid_1dx2d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { int nxgrid; nxgrid = create_xgrid_1dx2d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); return nxgrid; }; int create_xgrid_1dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { int nx1, ny1, nx2, ny2, nx1p, nx2p; int i1, j1, i2, j2, nxgrid, n; double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV]; double *area_in, *area_out, min_area; double *tmpx, *tmpy; nx1 = *nlon_in; ny1 = *nlat_in; nx2 = *nlon_out; ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; nx2p = nx2 + 1; area_in = (double *)malloc(nx1*ny1*sizeof(double)); area_out = (double *)malloc(nx2*ny2*sizeof(double)); tmpx = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double)); tmpy = (double *)malloc((nx1+1)*(ny1+1)*sizeof(double)); for(j1=0; j1<=ny1; j1++) for(i1=0; i1<=nx1; i1++) { tmpx[j1*nx1p+i1] = lon_in[i1]; tmpy[j1*nx1p+i1] = lat_in[j1]; } get_grid_area(nlon_in, nlat_in, tmpx, tmpy, area_in); get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out); free(tmpx); free(tmpy); for(j1=0; j1 MASK_THRESH ) { ll_lon = lon_in[i1]; ll_lat = lat_in[j1]; ur_lon = lon_in[i1+1]; ur_lat = lat_in[j1+1]; for(j2=0; j2=ur_lat) && (y_in[1]>=ur_lat) && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_out[j2*nx2p+i2]; x_in[1] = lon_out[j2*nx2p+i2+1]; x_in[2] = lon_out[(j2+1)*nx2p+i2+1]; x_in[3] = lon_out[(j2+1)*nx2p+i2]; n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); lon_in_avg = avgval_double(n_in, x_in); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if(xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = xarea; xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area_in); free(area_out); return nxgrid; }; /* create_xgrid_1dx2d_order2 */ /******************************************************************************* void create_xgrid_2dx1d_order1 This routine generate exchange grids between two grids for the first order conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds. mask is on grid lon_in/lat_in. *******************************************************************************/ int create_xgrid_2dx1d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) { int nxgrid; nxgrid = create_xgrid_2dx1d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, j_in, i_out, j_out, xgrid_area); return nxgrid; }; int create_xgrid_2dx1d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) { int nx1, ny1, nx2, ny2, nx1p, nx2p; int i1, j1, i2, j2, nxgrid; double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV]; double *area_in, *area_out, min_area; double *tmpx, *tmpy; nx1 = *nlon_in; ny1 = *nlat_in; nx2 = *nlon_out; ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; nx2p = nx2 + 1; area_in = (double *)malloc(nx1*ny1*sizeof(double)); area_out = (double *)malloc(nx2*ny2*sizeof(double)); tmpx = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double)); tmpy = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double)); for(j2=0; j2<=ny2; j2++) for(i2=0; i2<=nx2; i2++) { tmpx[j2*nx2p+i2] = lon_out[i2]; tmpy[j2*nx2p+i2] = lat_out[j2]; } get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in); get_grid_area(nlon_out, nlat_out, tmpx, tmpy, area_out); free(tmpx); free(tmpy); for(j2=0; j2 MASK_THRESH ) { int n_in, n_out; double Xarea; y_in[0] = lat_in[j1*nx1p+i1]; y_in[1] = lat_in[j1*nx1p+i1+1]; y_in[2] = lat_in[(j1+1)*nx1p+i1+1]; y_in[3] = lat_in[(j1+1)*nx1p+i1]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_in[j1*nx1p+i1]; x_in[1] = lon_in[j1*nx1p+i1+1]; x_in[2] = lon_in[(j1+1)*nx1p+i1+1]; x_in[3] = lon_in[(j1+1)*nx1p+i1]; n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { Xarea = poly_area ( x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if( Xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = Xarea; i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area_in); free(area_out); return nxgrid; }; /* create_xgrid_2dx1d_order1 */ /******************************************************************************** void create_xgrid_2dx1d_order2 This routine generate exchange grids between two grids for the second order conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_out,lat_out are 1-D grid bounds, lon_in,lat_in are geographic grid location of grid cell bounds. mask is on grid lon_in/lat_in. ********************************************************************************/ int create_xgrid_2dx1d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { int nxgrid; nxgrid = create_xgrid_2dx1d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); return nxgrid; }; int create_xgrid_2dx1d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { int nx1, ny1, nx2, ny2, nx1p, nx2p; int i1, j1, i2, j2, nxgrid, n; double ll_lon, ll_lat, ur_lon, ur_lat, x_in[MV], y_in[MV], x_out[MV], y_out[MV]; double *tmpx, *tmpy; double *area_in, *area_out, min_area; double lon_in_avg; nx1 = *nlon_in; ny1 = *nlat_in; nx2 = *nlon_out; ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; nx2p = nx2 + 1; area_in = (double *)malloc(nx1*ny1*sizeof(double)); area_out = (double *)malloc(nx2*ny2*sizeof(double)); tmpx = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double)); tmpy = (double *)malloc((nx2+1)*(ny2+1)*sizeof(double)); for(j2=0; j2<=ny2; j2++) for(i2=0; i2<=nx2; i2++) { tmpx[j2*nx2p+i2] = lon_out[i2]; tmpy[j2*nx2p+i2] = lat_out[j2]; } get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in); get_grid_area(nlon_out, nlat_out, tmpx, tmpy, area_out); free(tmpx); free(tmpy); for(j2=0; j2 MASK_THRESH ) { int n_in, n_out; double xarea; y_in[0] = lat_in[j1*nx1p+i1]; y_in[1] = lat_in[j1*nx1p+i1+1]; y_in[2] = lat_in[(j1+1)*nx1p+i1+1]; y_in[3] = lat_in[(j1+1)*nx1p+i1]; if ( (y_in[0]<=ll_lat) && (y_in[1]<=ll_lat) && (y_in[2]<=ll_lat) && (y_in[3]<=ll_lat) ) continue; if ( (y_in[0]>=ur_lat) && (y_in[1]>=ur_lat) && (y_in[2]>=ur_lat) && (y_in[3]>=ur_lat) ) continue; x_in[0] = lon_in[j1*nx1p+i1]; x_in[1] = lon_in[j1*nx1p+i1+1]; x_in[2] = lon_in[(j1+1)*nx1p+i1+1]; x_in[3] = lon_in[(j1+1)*nx1p+i1]; n_in = fix_lon(x_in, y_in, 4, (ll_lon+ur_lon)/2); lon_in_avg = avgval_double(n_in, x_in); if ( (n_out = clip ( x_in, y_in, n_in, ll_lon, ll_lat, ur_lon, ur_lat, x_out, y_out )) > 0 ) { xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if(xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = xarea; xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area_in); free(area_out); return nxgrid; }; /* create_xgrid_2dx1d_order2 */ /******************************************************************************* void create_xgrid_2DX2D_order1 This routine generate exchange grids between two grids for the first order conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds. mask is on grid lon_in/lat_in. *******************************************************************************/ #ifndef __AIX int create_xgrid_2dx2d_order1_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) { int nxgrid; nxgrid = create_xgrid_2dx2d_order1(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, j_in, i_out, j_out, xgrid_area); return nxgrid; }; #endif int create_xgrid_2dx2d_order1(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area) { int nx1, ny1, nx2, ny2, nx1p, nx2p, ny1p, ny2p, nxgrid, n1_in, n2_in; int n0, n1, n2, n3, i1, j1, i2, j2, l; double x1_in[MV], y1_in[MV], x2_in[MV], y2_in[MV], x_out[MV], y_out[MV]; double lon_in_min, lon_out_min, lon_in_max, lon_out_max, lat_in_min, lat_in_max, lat_out_min, lat_out_max; double lon_in_avg; double *area_in, *area_out, min_area; nx1 = *nlon_in; ny1 = *nlat_in; nx2 = *nlon_out; ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; nx2p = nx2 + 1; area_in = (double *)malloc(nx1*ny1*sizeof(double)); area_out = (double *)malloc(nx2*ny2*sizeof(double)); get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in); get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out); for(j1=0; j1 MASK_THRESH ) { n0 = j1*nx1p+i1; n1 = j1*nx1p+i1+1; n2 = (j1+1)*nx1p+i1+1; n3 = (j1+1)*nx1p+i1; x1_in[0] = lon_in[n0]; y1_in[0] = lat_in[n0]; x1_in[1] = lon_in[n1]; y1_in[1] = lat_in[n1]; x1_in[2] = lon_in[n2]; y1_in[2] = lat_in[n2]; x1_in[3] = lon_in[n3]; y1_in[3] = lat_in[n3]; lat_in_min = minval_double(4, y1_in); lat_in_max = maxval_double(4, y1_in); n1_in = fix_lon(x1_in, y1_in, 4, M_PI); lon_in_min = minval_double(n1_in, x1_in); lon_in_max = maxval_double(n1_in, x1_in); lon_in_avg = avgval_double(n1_in, x1_in); for(j2=0; j2= lat_in_max || lat_out_max <= lat_in_min ) continue; n2_in = fix_lon(x2_in, y2_in, 4, lon_in_avg); lon_out_min = minval_double(n2_in, x2_in); lon_out_max = maxval_double(n2_in, x2_in); /* x2_in should in the same range as lon_in_in after lon_fix, so no need to consider cyclic condition */ if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue; if ( (n_out = clip_2dx2d ( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) { Xarea = poly_area (x_out, y_out, n_out) * mask_in[j1*nx1+i1]; min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if( Xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = Xarea; i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area_in); free(area_out); return nxgrid; };/* get_xgrid_2Dx2D_order1 */ /******************************************************************************** void create_xgrid_2dx1d_order2 This routine generate exchange grids between two grids for the second order conservative interpolation. nlon_in,nlat_in,nlon_out,nlat_out are the size of the grid cell and lon_in,lat_in, lon_out,lat_out are geographic grid location of grid cell bounds. mask is on grid lon_in/lat_in. ********************************************************************************/ #ifndef __AIX int create_xgrid_2dx2d_order2_(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { int nxgrid; nxgrid = create_xgrid_2dx2d_order2(nlon_in, nlat_in, nlon_out, nlat_out, lon_in, lat_in, lon_out, lat_out, mask_in, i_in, j_in, i_out, j_out, xgrid_area, xgrid_clon, xgrid_clat); return nxgrid; }; #endif int create_xgrid_2dx2d_order2(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out, const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out, const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, double *xgrid_area, double *xgrid_clon, double *xgrid_clat) { int nx1, nx2, ny1, ny2, nx1p, nx2p, ny1p, ny2p, nxgrid, n1_in, n2_in; int n0, n1, n2, n3, i1, j1, i2, j2, l, n; double x1_in[MV], y1_in[MV], x2_in[MV], y2_in[MV], x_out[MV], y_out[MV]; double lon_in_min, lon_out_min, lon_in_max, lon_out_max, lat_in_min, lat_in_max, lat_out_min, lat_out_max; double lon_in_avg, xctrlon, xctrlat; double *area_in, *area_out, min_area; nx1 = *nlon_in; ny1 = *nlat_in; nx2 = *nlon_out; ny2 = *nlat_out; nxgrid = 0; nx1p = nx1 + 1; nx2p = nx2 + 1; ny1p = ny1 + 1; ny2p = ny2 + 1; area_in = (double *)malloc(nx1*ny1*sizeof(double)); area_out = (double *)malloc(nx2*ny2*sizeof(double)); get_grid_area(nlon_in, nlat_in, lon_in, lat_in, area_in); get_grid_area(nlon_out, nlat_out, lon_out, lat_out, area_out); for(j1=0; j1 MASK_THRESH ) { n0 = j1*nx1p+i1; n1 = j1*nx1p+i1+1; n2 = (j1+1)*nx1p+i1+1; n3 = (j1+1)*nx1p+i1; x1_in[0] = lon_in[n0]; y1_in[0] = lat_in[n0]; x1_in[1] = lon_in[n1]; y1_in[1] = lat_in[n1]; x1_in[2] = lon_in[n2]; y1_in[2] = lat_in[n2]; x1_in[3] = lon_in[n3]; y1_in[3] = lat_in[n3]; lat_in_min = minval_double(4, y1_in); lat_in_max = maxval_double(4, y1_in); n1_in = fix_lon(x1_in, y1_in, 4, M_PI); lon_in_min = minval_double(n1_in, x1_in); lon_in_max = maxval_double(n1_in, x1_in); lon_in_avg = avgval_double(n1_in, x1_in); for(j2=0; j2= lat_in_max || lat_out_max <= lat_in_min ) continue; n2_in = fix_lon(x2_in, y2_in, 4, lon_in_avg); lon_out_min = minval_double(n2_in, x2_in); lon_out_max = maxval_double(n2_in, x2_in); /* x2_in should in the same range as x1_in after lon_fix, so no need to consider cyclic condition */ if(lon_out_min >= lon_in_max || lon_out_max <= lon_in_min ) continue; if ( (n_out = clip_2dx2d( x1_in, y1_in, n1_in, x2_in, y2_in, n2_in, x_out, y_out )) > 0) { xarea = poly_area (x_out, y_out, n_out ) * mask_in[j1*nx1+i1]; min_area = min(area_in[j1*nx1+i1], area_out[j2*nx2+i2]); if( xarea/min_area > AREA_RATIO_THRESH ) { xgrid_area[nxgrid] = xarea; xgrid_clon[nxgrid] = poly_ctrlon(x_out, y_out, n_out, lon_in_avg); xgrid_clat[nxgrid] = poly_ctrlat (x_out, y_out, n_out ); i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area_in); free(area_out); return nxgrid; };/* get_xgrid_2Dx2D_order2 */ /******************************************************************************* Sutherland-Hodgeman algorithm sequentially clips parts outside 4 boundaries *******************************************************************************/ int clip(const double lon_in[], const double lat_in[], int n_in, double ll_lon, double ll_lat, double ur_lon, double ur_lat, double lon_out[], double lat_out[]) { double x_tmp[MV], y_tmp[MV], x_last, y_last; int i_in, i_out, n_out, inside_last, inside; /* clip polygon with LEFT boundary - clip V_IN to V_TMP */ x_last = lon_in[n_in-1]; y_last = lat_in[n_in-1]; inside_last = (x_last >= ll_lon); for (i_in=0,i_out=0;i_in= ll_lon))!=inside_last) { x_tmp[i_out] = ll_lon; y_tmp[i_out++] = y_last + (ll_lon - x_last) * (lat_in[i_in] - y_last) / (lon_in[i_in] - x_last); } /* if "to" point is right of LEFT boundary, output it */ if (inside) { x_tmp[i_out] = lon_in[i_in]; y_tmp[i_out++] = lat_in[i_in]; } x_last = lon_in[i_in]; y_last = lat_in[i_in]; inside_last = inside; } if (!(n_out=i_out)) return(0); /* clip polygon with RIGHT boundary - clip V_TMP to V_OUT */ x_last = x_tmp[n_out-1]; y_last = y_tmp[n_out-1]; inside_last = (x_last <= ur_lon); for (i_in=0,i_out=0;i_in= ll_lat); for (i_in=0,i_out=0;i_in= ll_lat))!=inside_last) { y_tmp[i_out] = ll_lat; x_tmp[i_out++] = x_last + (ll_lat - y_last) * (lon_out[i_in] - x_last) / (lat_out[i_in] - y_last); } /* if "to" point is above BOTTOM boundary, output it */ if (inside) { x_tmp[i_out] = lon_out[i_in]; y_tmp[i_out++] = lat_out[i_in]; } x_last = lon_out[i_in]; y_last = lat_out[i_in]; inside_last = inside; } if (!(n_out=i_out)) return(0); /* clip polygon with TOP boundary - clip V_TMP to V_OUT */ x_last = x_tmp[n_out-1]; y_last = y_tmp[n_out-1]; inside_last = (y_last <= ur_lat); for (i_in=0,i_out=0;i_in and should not parallel to the line between and may need to consider truncation error */ dy1 = y1_1-y1_0; dy2 = y2_1-y2_0; dx1 = x1_1-x1_0; dx2 = x2_1-x2_0; ds1 = y1_0*x1_1 - y1_1*x1_0; ds2 = y2_0*x2_1 - y2_1*x2_0; determ = dy2*dx1 - dy1*dx2; if(fabs(determ) < EPSLN30) { error_handler("the line between and should not parallel to " "the line between and "); } lon_out[i_out] = (dx2*ds1 - dx1*ds2)/determ; lat_out[i_out++] = (dy2*ds1 - dy1*ds2)/determ; } if(inside) { lon_out[i_out] = x1_1; lat_out[i_out++] = y1_1; } x1_0 = x1_1; y1_0 = y1_1; inside_last = inside; } if(!(n_out=i_out)) return 0; for(i1=0; i1 MASK_THRESH ) { /* clockwise */ n0 = j1*nx1p+i1; n1 = (j1+1)*nx1p+i1; n2 = (j1+1)*nx1p+i1+1; n3 = j1*nx1p+i1+1; x1_in[0] = x1[n0]; y1_in[0] = y1[n0]; z1_in[0] = z1[n0]; x1_in[1] = x1[n1]; y1_in[1] = y1[n1]; z1_in[1] = z1[n1]; x1_in[2] = x1[n2]; y1_in[2] = y1[n2]; z1_in[2] = z1[n2]; x1_in[3] = x1[n3]; y1_in[3] = y1[n3]; z1_in[3] = z1[n3]; for(j2=0; j2 0) { xarea = great_circle_area ( n_out, x_out, y_out, z_out ) * mask_in[j1*nx1+i1]; min_area = min(area1[j1*nx1+i1], area2[j2*nx2+i2]); if( xarea/min_area > AREA_RATIO_THRESH ) { #ifdef debug_test_create_xgrid printf("(i2,j2)=(%d,%d), (i1,j1)=(%d,%d), xarea=%g\n", i2, j2, i1, j1, xarea); #endif xgrid_area[nxgrid] = xarea; xgrid_clon[nxgrid] = 0; /*z1l: will be developed very soon */ xgrid_clat[nxgrid] = 0; i_in[nxgrid] = i1; j_in[nxgrid] = j1; i_out[nxgrid] = i2; j_out[nxgrid] = j2; ++nxgrid; if(nxgrid > MAXXGRID) error_handler("nxgrid is greater than MAXXGRID, increase MAXXGRID"); } } } } free(area1); free(area2); free(x1); free(y1); free(z1); free(x2); free(y2); free(z2); return nxgrid; };/* create_xgrid_great_circle */ /******************************************************************************* Revise Sutherland-Hodgeman algorithm to find the vertices of the overlapping between any two grid boxes. It return the number of vertices for the exchange grid. Each edge of grid box is a part of great circle. All the points are cartesian coordinates. Here we are assuming each polygon is convex. RANGE_CHECK_CRITERIA is used to determine if the two grid boxes are possible to be overlap. The size should be between 0 and 0.5. The larger the range_check_criteria, the more expensive of the computatioin. When the value is close to 0, some small exchange grid might be lost. Suggest to use value 0.05 for C48. *******************************************************************************/ int clip_2dx2d_great_circle(const double x1_in[], const double y1_in[], const double z1_in[], int n1_in, const double x2_in[], const double y2_in[], const double z2_in [], int n2_in, double x_out[], double y_out[], double z_out[]) { struct Node *subjList=NULL; struct Node *clipList=NULL; struct Node *grid1List=NULL; struct Node *grid2List=NULL; struct Node *intersectList=NULL; struct Node *polyList=NULL; struct Node *curList=NULL; struct Node *firstIntersect=NULL, *curIntersect=NULL; struct Node *temp1=NULL, *temp2=NULL, *temp=NULL; int i1, i2, i1p, i2p, i2p2, npts1, npts2; int nintersect, n_out; int maxiter1, maxiter2, iter1, iter2; int found1, found2, curListNum; int has_inbound, inbound; double pt1[MV][3], pt2[MV][3]; double *p1_0=NULL, *p1_1=NULL; double *p2_0=NULL, *p2_1=NULL, *p2_2=NULL; double intersect[3]; double u1, u2; double min_x1, max_x1, min_y1, max_y1, min_z1, max_z1; double min_x2, max_x2, min_y2, max_y2, min_z2, max_z2; static int first_call=1; /* first check the min and max of (x1_in, y1_in, z1_in) with (x2_in, y2_in, z2_in) */ min_x1 = minval_double(n1_in, x1_in); max_x2 = maxval_double(n2_in, x2_in); if(min_x1 >= max_x2+RANGE_CHECK_CRITERIA) return 0; max_x1 = maxval_double(n1_in, x1_in); min_x2 = minval_double(n2_in, x2_in); if(min_x2 >= max_x1+RANGE_CHECK_CRITERIA) return 0; min_y1 = minval_double(n1_in, y1_in); max_y2 = maxval_double(n2_in, y2_in); if(min_y1 >= max_y2+RANGE_CHECK_CRITERIA) return 0; max_y1 = maxval_double(n1_in, y1_in); min_y2 = minval_double(n2_in, y2_in); if(min_y2 >= max_y1+RANGE_CHECK_CRITERIA) return 0; min_z1 = minval_double(n1_in, z1_in); max_z2 = maxval_double(n2_in, z2_in); if(min_z1 >= max_z2+RANGE_CHECK_CRITERIA) return 0; max_z1 = maxval_double(n1_in, z1_in); min_z2 = minval_double(n2_in, z2_in); if(min_z2 >= max_z1+RANGE_CHECK_CRITERIA) return 0; rewindList(); grid1List = getNext(); grid2List = getNext(); intersectList = getNext(); polyList = getNext(); /* insert points into SubjList and ClipList */ for(i1=0; i1isInside = 1; else temp->isInside = 0; temp = getNextNode(temp); } #ifdef debug_test_create_xgrid printf("\nNOTE from clip_2dx2d_great_circle: begin to set inside value of grid2List\n"); #endif /* check if grid2List is inside grid1List */ temp = grid2List; while(temp) { if(insidePolygon(temp, grid1List)) temp->isInside = 1; else temp->isInside = 0; temp = getNextNode(temp); } /* make sure the grid box is clockwise */ /*make sure each polygon is convex, which is equivalent that the great_circle_area is positive */ if( gridArea(grid1List) <= 0 ) error_handler("create_xgrid.c(clip_2dx2d_great_circle): grid box 1 is not convex"); if( gridArea(grid2List) <= 0 ) error_handler("create_xgrid.c(clip_2dx2d_great_circle): grid box 2 is not convex"); #ifdef debug_test_create_xgrid printNode(grid1List, "grid1List"); printNode(grid2List, "grid2List"); #endif /* get the coordinates from grid1List and grid2List. Please not npts1 might not equal n1_in, npts2 might not equal n2_in because of pole */ temp = grid1List; for(i1=0; i1Next; } temp = grid2List; for(i2=0; i2Next; } firstIntersect=getNext(); curIntersect = getNext(); #ifdef debug_test_create_xgrid printf("\n\n************************ Start line_intersect_2D_3D ******************************\n"); #endif /* first find all the intersection points */ nintersect = 0; for(i1=0; i1 1) { getFirstInbound(intersectList, firstIntersect); if(firstIntersect->initialized) { has_inbound = 1; } } /* when has_inbound == 0, get the grid1List and grid2List */ if( !has_inbound && nintersect > 1) { setInbound(intersectList, grid1List); getFirstInbound(intersectList, firstIntersect); if(firstIntersect->initialized) has_inbound = 1; } /* if has_inbound = 1, find the overlapping */ n_out = 0; if(has_inbound) { maxiter1 = nintersect; #ifdef debug_test_create_xgrid printf("\nNOTE from clip_2dx2d_great_circle: number of intersect is %d\n", nintersect); printf("\n size of grid2List is %d, size of grid1List is %d\n", length(grid2List), length(grid1List)); printNode(intersectList, "beginning intersection list"); printNode(grid2List, "beginning clip list"); printNode(grid1List, "beginning subj list"); printf("\n************************ End line_intersect_2D_3D **********************************\n\n"); #endif temp1 = getNode(grid1List, *firstIntersect); if( temp1 == NULL) { double lon[10], lat[10]; int i; xyz2latlon(n1_in, x1_in, y1_in, z1_in, lon, lat); for(i=0; i< n1_in; i++) printf("lon1 = %g, lat1 = %g\n", lon[i]*R2D, lat[i]*R2D); printf("\n"); xyz2latlon(n2_in, x2_in, y2_in, z2_in, lon, lat); for(i=0; i< n2_in; i++) printf("lon2 = %g, lat2 = %g\n", lon[i]*R2D, lat[i]*R2D); printf("\n"); error_handler("firstIntersect is not in the grid1List"); } addNode(polyList, *firstIntersect); nintersect--; #ifdef debug_test_create_xgrid printNode(polyList, "polyList at stage 1"); #endif /* Loop over the grid1List and grid2List to find again the firstIntersect */ curList = grid1List; curListNum = 0; /* Loop through curList to find the next intersection, the loop will end when come back to firstIntersect */ copyNode(curIntersect, *firstIntersect); iter1 = 0; found1 = 0; while( iter1 < maxiter1 ) { #ifdef debug_test_create_xgrid printf("\n----------- At iteration = %d\n\n", iter1+1 ); printNode(curIntersect, "curIntersect at the begining of iter1"); #endif /* find the curIntersect in curList and get the next intersection points */ temp1 = getNode(curList, *curIntersect); temp2 = temp1->Next; if( temp2 == NULL ) temp2 = curList; maxiter2 = length(curList); found2 = 0; iter2 = 0; /* Loop until find the next intersection */ while( iter2 < maxiter2 ) { int temp2IsIntersect; temp2IsIntersect = 0; if( isIntersect( *temp2 ) ) { /* copy the point and switch to the grid2List */ struct Node *temp3; /* first check if temp2 is the firstIntersect */ if( sameNode( *temp2, *firstIntersect) ) { found1 = 1; break; } temp3 = temp2->Next; if( temp3 == NULL) temp3 = curList; if( temp3 == NULL) error_handler("creat_xgrid.c: temp3 can not be NULL"); found2 = 1; /* if next node is inside or an intersection, need to keep on curList */ temp2IsIntersect = 1; if( isIntersect(*temp3) || (temp3->isInside == 1) ) found2 = 0; } if(found2) { copyNode(curIntersect, *temp2); break; } else { addNode(polyList, *temp2); #ifdef debug_test_create_xgrid printNode(polyList, "polyList at stage 2"); #endif if(temp2IsIntersect) { nintersect--; } } temp2 = temp2->Next; if( temp2 == NULL ) temp2 = curList; iter2 ++; } if(found1) break; if( !found2 ) error_handler(" not found the next intersection "); /* if find the first intersection, the poly found */ if( sameNode( *curIntersect, *firstIntersect) ) { found1 = 1; break; } /* add curIntersect to polyList and remove it from intersectList and curList */ addNode(polyList, *curIntersect); #ifdef debug_test_create_xgrid printNode(polyList, "polyList at stage 3"); #endif nintersect--; /* switch curList */ if( curListNum == 0) { curList = grid2List; curListNum = 1; } else { curList = grid1List; curListNum = 0; } iter1++; } if(!found1) error_handler("not return back to the first intersection"); /* currently we are only clipping convex polygon to convex polygon */ if( nintersect > 0) error_handler("After clipping, nintersect should be 0"); /* copy the polygon to x_out, y_out, z_out */ temp1 = polyList; while (temp1 != NULL) { getCoordinate(*temp1, x_out+n_out, y_out+n_out, z_out+n_out); temp1 = temp1->Next; n_out++; } /* if(n_out < 3) error_handler(" The clipped region has < 3 vertices"); */ if( n_out < 3) n_out = 0; #ifdef debug_test_create_xgrid printNode(polyList, "polyList after clipping"); #endif } /* check if grid1 is inside grid2 */ if(n_out==0){ /* first check number of points in grid1 is inside grid2 */ int n, n1in2; /* One possible is that grid1List is inside grid2List */ #ifdef debug_test_create_xgrid printf("\nNOTE from clip_2dx2d_great_circle: check if grid1 is inside grid2\n"); #endif n1in2 = 0; temp = grid1List; while(temp) { if(temp->intersect != 1) { #ifdef debug_test_create_xgrid printf("grid1->isInside = %d\n", temp->isInside); #endif if( temp->isInside == 1) n1in2++; } temp = getNextNode(temp); } if(npts1==n1in2) { /* grid1 is inside grid2 */ n_out = npts1; n = 0; temp = grid1List; while( temp ) { getCoordinate(*temp, &x_out[n], &y_out[n], &z_out[n]); n++; temp = getNextNode(temp); } } if(n_out>0) return n_out; } /* check if grid2List is inside grid1List */ if(n_out ==0){ int n, n2in1; #ifdef debug_test_create_xgrid printf("\nNOTE from clip_2dx2d_great_circle: check if grid2 is inside grid1\n"); #endif temp = grid2List; n2in1 = 0; while(temp) { if(temp->intersect != 1) { #ifdef debug_test_create_xgrid printf("grid2->isInside = %d\n", temp->isInside); #endif if( temp->isInside == 1) n2in1++; } temp = getNextNode(temp); } if(npts2==n2in1) { /* grid2 is inside grid1 */ n_out = npts2; n = 0; temp = grid2List; while( temp ) { getCoordinate(*temp, &x_out[n], &y_out[n], &z_out[n]); n++; temp = getNextNode(temp); } } } return n_out; } /* Intersects between the line a and the seqment s where both line and segment are great circle lines on the sphere represented by 3D cartesian points. [sin sout] are the ends of a line segment returns true if the lines could be intersected, false otherwise. inbound means the direction of (a1,a2) go inside or outside of (q1,q2,q3) */ int line_intersect_2D_3D(double *a1, double *a2, double *q1, double *q2, double *q3, double *intersect, double *u_a, double *u_q, int *inbound){ /* Do this intersection by reprsenting the line a1 to a2 as a plane through the two line points and the origin of the sphere (0,0,0). This is the definition of a great circle arc. */ double plane[9]; double plane_p[2]; double u; double p1[3], v1[3], v2[3]; double c1[3], c2[3], c3[3]; double coincident, sense, norm; int i; int is_inter1, is_inter2; *inbound = 0; /* first check if any vertices are the same */ if(samePoint(a1[0], a1[1], a1[2], q1[0], q1[1], q1[2])) { *u_a = 0; *u_q = 0; intersect[0] = a1[0]; intersect[1] = a1[1]; intersect[2] = a1[2]; #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound); #endif return 1; } else if (samePoint(a1[0], a1[1], a1[2], q2[0], q2[1], q2[2])) { *u_a = 0; *u_q = 1; intersect[0] = a1[0]; intersect[1] = a1[1]; intersect[2] = a1[2]; #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound); #endif return 1; } else if(samePoint(a2[0], a2[1], a2[2], q1[0], q1[1], q1[2])) { #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound); #endif *u_a = 1; *u_q = 0; intersect[0] = a2[0]; intersect[1] = a2[1]; intersect[2] = a2[2]; return 1; } else if (samePoint(a2[0], a2[1], a2[2], q2[0], q2[1], q2[2])) { #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f, u_q=%19.15f, inbound=%d\n", *u_a, *u_q, *inbound); #endif *u_a = 1; *u_q = 1; intersect[0] = a2[0]; intersect[1] = a2[1]; intersect[2] = a2[2]; return 1; } /* Load points defining plane into variable (these are supposed to be in counterclockwise order) */ plane[0]=q1[0]; plane[1]=q1[1]; plane[2]=q1[2]; plane[3]=q2[0]; plane[4]=q2[1]; plane[5]=q2[2]; plane[6]=0.0; plane[7]=0.0; plane[8]=0.0; /* Intersect the segment with the plane */ is_inter1 = intersect_tri_with_line(plane, a1, a2, plane_p, u_a); if(!is_inter1) return 0; if(fabs(*u_a) < EPSLN8) *u_a = 0; if(fabs(*u_a-1) < EPSLN8) *u_a = 1; #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: u_a = %19.15f\n", *u_a); #endif if( (*u_a < 0) || (*u_a > 1) ) return 0; /* Load points defining plane into variable (these are supposed to be in counterclockwise order) */ plane[0]=a1[0]; plane[1]=a1[1]; plane[2]=a1[2]; plane[3]=a2[0]; plane[4]=a2[1]; plane[5]=a2[2]; plane[6]=0.0; plane[7]=0.0; plane[8]=0.0; /* Intersect the segment with the plane */ is_inter2 = intersect_tri_with_line(plane, q1, q2, plane_p, u_q); if(!is_inter2) return 0; if(fabs(*u_q) < EPSLN8) *u_q = 0; if(fabs(*u_q-1) < EPSLN8) *u_q = 1; #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: u_q = %19.15f\n", *u_q); #endif if( (*u_q < 0) || (*u_q > 1) ) return 0; u =*u_a; /* The two planes are coincidental */ vect_cross(a1, a2, c1); vect_cross(q1, q2, c2); vect_cross(c1, c2, c3); coincident = metric(c3); if(fabs(coincident) < EPSLN30) return 0; /* Calculate point of intersection */ intersect[0]=a1[0] + u*(a2[0]-a1[0]); intersect[1]=a1[1] + u*(a2[1]-a1[1]); intersect[2]=a1[2] + u*(a2[2]-a1[2]); norm = metric( intersect ); for(i = 0; i < 3; i ++) intersect[i] /= norm; /* when u_q =0 or u_q =1, the following could not decide the inbound value */ if(*u_q != 0 && *u_q != 1){ p1[0] = a2[0]-a1[0]; p1[1] = a2[1]-a1[1]; p1[2] = a2[2]-a1[2]; v1[0] = q2[0]-q1[0]; v1[1] = q2[1]-q1[1]; v1[2] = q2[2]-q1[2]; v2[0] = q3[0]-q2[0]; v2[1] = q3[1]-q2[1]; v2[2] = q3[2]-q2[2]; vect_cross(v1, v2, c1); vect_cross(v1, p1, c2); sense = dot(c1, c2); *inbound = 1; if(sense > 0) *inbound = 2; /* v1 going into v2 in CCW sense */ } #ifdef debug_test_create_xgrid printf("\nNOTE from line_intersect_2D_3D: inbound=%d\n", *inbound); #endif return 1; } /*------------------------------------------------------------------------------ double poly_ctrlat(const double x[], const double y[], int n) This routine is used to calculate the latitude of the centroid ---------------------------------------------------------------------------*/ double poly_ctrlat(const double x[], const double y[], int n) { double ctrlat = 0.0; int i; for (i=0;i M_PI) dx = dx - 2.0*M_PI; if(dx < -M_PI) dx = dx + 2.0*M_PI; if ( fabs(hdy)< SMALL_VALUE ) /* cheap area calculation along latitude */ ctrlat -= dx*(2*cos(avg_y) + lat2*sin(avg_y) - cos(lat1) ); else ctrlat -= dx*( (sin(hdy)/hdy)*(2*cos(avg_y) + lat2*sin(avg_y)) - cos(lat1) ); } return (ctrlat*RADIUS*RADIUS); }; /* poly_ctrlat */ /*------------------------------------------------------------------------------ double poly_ctrlon(const double x[], const double y[], int n, double clon) This routine is used to calculate the lontitude of the centroid. ---------------------------------------------------------------------------*/ double poly_ctrlon(const double x[], const double y[], int n, double clon) { double ctrlon = 0.0; int i; for (i=0;i M_PI) dphi = dphi - 2.0*M_PI; if(dphi < -M_PI) dphi = dphi + 2.0*M_PI; dphi1 = phi1 - clon; if( dphi1 > M_PI) dphi1 -= 2.0*M_PI; if( dphi1 <-M_PI) dphi1 += 2.0*M_PI; dphi2 = phi2 -clon; if( dphi2 > M_PI) dphi2 -= 2.0*M_PI; if( dphi2 <-M_PI) dphi2 += 2.0*M_PI; if(fabs(dphi2 -dphi1) < M_PI) { ctrlon -= dphi * (dphi1*f1+dphi2*f2)/2.0; } else { if(dphi1 > 0.0) fac = M_PI; else fac = -M_PI; fint = f1 + (f2-f1)*(fac-dphi1)/fabs(dphi); ctrlon -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2 + 0.5*fac*(dphi1+dphi2)*fint; } } return (ctrlon*RADIUS*RADIUS); }; /* poly_ctrlon */ /* ----------------------------------------------------------------------------- double box_ctrlat(double ll_lon, double ll_lat, double ur_lon, double ur_lat) This routine is used to calculate the latitude of the centroid. ---------------------------------------------------------------------------*/ double box_ctrlat(double ll_lon, double ll_lat, double ur_lon, double ur_lat) { double dphi = ur_lon-ll_lon; double ctrlat; if(dphi > M_PI) dphi = dphi - 2.0*M_PI; if(dphi < -M_PI) dphi = dphi + 2.0*M_PI; ctrlat = dphi*(cos(ur_lat) + ur_lat*sin(ur_lat)-(cos(ll_lat) + ll_lat*sin(ll_lat))); return (ctrlat*RADIUS*RADIUS); }; /* box_ctrlat */ /*------------------------------------------------------------------------------ double box_ctrlon(double ll_lon, double ll_lat, double ur_lon, double ur_lat, double clon) This routine is used to calculate the lontitude of the centroid ----------------------------------------------------------------------------*/ double box_ctrlon(double ll_lon, double ll_lat, double ur_lon, double ur_lat, double clon) { double phi1, phi2, dphi, lat1, lat2, dphi1, dphi2; double f1, f2, fac, fint; double ctrlon = 0.0; int i; for( i =0; i<2; i++) { if(i == 0) { phi1 = ur_lon; phi2 = ll_lon; lat1 = lat2 = ll_lat; } else { phi1 = ll_lon; phi2 = ur_lon; lat1 = lat2 = ur_lat; } dphi = phi1 - phi2; f1 = 0.5*(cos(lat1)*sin(lat1)+lat1); f2 = 0.5*(cos(lat2)*sin(lat2)+lat2); if(dphi > M_PI) dphi = dphi - 2.0*M_PI; if(dphi < -M_PI) dphi = dphi + 2.0*M_PI; /* make sure the center is in the same grid box. */ dphi1 = phi1 - clon; if( dphi1 > M_PI) dphi1 -= 2.0*M_PI; if( dphi1 <-M_PI) dphi1 += 2.0*M_PI; dphi2 = phi2 -clon; if( dphi2 > M_PI) dphi2 -= 2.0*M_PI; if( dphi2 <-M_PI) dphi2 += 2.0*M_PI; if(fabs(dphi2 -dphi1) < M_PI) { ctrlon -= dphi * (dphi1*f1+dphi2*f2)/2.0; } else { if(dphi1 > 0.0) fac = M_PI; else fac = -M_PI; fint = f1 + (f2-f1)*(fac-dphi1)/fabs(dphi); ctrlon -= 0.5*dphi1*(dphi1-fac)*f1 - 0.5*dphi2*(dphi2+fac)*f2 + 0.5*fac*(dphi1+dphi2)*fint; } } return (ctrlon*RADIUS*RADIUS); } /* box_ctrlon */ /******************************************************************************* double grid_box_radius(double *x, double *y, double *z, int n); Find the radius of the grid box, the radius is defined the maximum distance between any two vertices *******************************************************************************/ double grid_box_radius(const double *x, const double *y, const double *z, int n) { double radius; int i, j; radius = 0; for(i=0; i is the outward edge normal from vertex to . is the vector from to . if Inner produce * > 0, outside, otherwise inside. inner product value = 0 also treate as inside. *******************************************************************************/ int inside_edge(double x0, double y0, double x1, double y1, double x, double y) { const double SMALL = 1.e-12; double product; product = ( x-x0 )*(y1-y0) + (x0-x1)*(y-y0); return (product<=SMALL) ? 1:0; }; /* inside_edge */ /* The following is a test program to test subroutines in create_xgrid.c */ #ifdef test_create_xgrid #include "create_xgrid.h" #include #define D2R (M_PI/180) #define R2D (180/M_PI) #define MAXPOINT 1000 int main(int argc, char* argv[]) { double lon1_in[MAXPOINT], lat1_in[MAXPOINT]; double lon2_in[MAXPOINT], lat2_in[MAXPOINT]; double x1_in[MAXPOINT], y1_in[MAXPOINT], z1_in[MAXPOINT]; double x2_in[MAXPOINT], y2_in[MAXPOINT], z2_in[MAXPOINT]; double lon_out[20], lat_out[20]; double x_out[20], y_out[20], z_out[20]; int n1_in, n2_in, n_out, i, j; int nlon1=0, nlat1=0, nlon2=0, nlat2=0; int n; int ntest = 11; for(n=11; n<=ntest; n++) { switch (n) { case 1: /**************************************************************** test clip_2dx2d_great_cirle case 1: box 1: (20,10), (20,12), (22,12), (22,10) box 2: (21,11), (21,14), (24,14), (24,11) out : (21, 12.0018), (22, 12), (22, 11.0033), (21, 11) ****************************************************************/ n1_in = 4; n2_in = 4; /* first a simple lat-lon grid box to clip another lat-lon grid box */ lon1_in[0] = 20; lat1_in[0] = 10; lon1_in[1] = 20; lat1_in[1] = 12; lon1_in[2] = 22; lat1_in[2] = 12; lon1_in[3] = 22; lat1_in[3] = 10; lon2_in[0] = 21; lat2_in[0] = 11; lon2_in[1] = 21; lat2_in[1] = 14; lon2_in[2] = 24; lat2_in[2] = 14; lon2_in[3] = 24; lat2_in[3] = 11; break; case 2: /**************************************************************** test clip_2dx2d_great_cirle case 2: two identical box box 1: (20,10), (20,12), (22,12), (22,10) box 2: (20,10), (20,12), (22,12), (22,10) out : (20,10), (20,12), (22,12), (22,10) ****************************************************************/ lon1_in[0] = 20; lat1_in[0] = 10; lon1_in[1] = 20; lat1_in[1] = 12; lon1_in[2] = 22; lat1_in[2] = 12; lon1_in[3] = 22; lat1_in[3] = 10; for(i=0; i 10 ) { int nxgrid; int *i1, *j1, *i2, *j2; double *xarea, *xclon, *xclat, *mask1; mask1 = (double *)malloc(nlon1*nlat1*sizeof(double)); i1 = (int *)malloc(MAXXGRID*sizeof(int)); j1 = (int *)malloc(MAXXGRID*sizeof(int)); i2 = (int *)malloc(MAXXGRID*sizeof(int)); j2 = (int *)malloc(MAXXGRID*sizeof(int)); xarea = (double *)malloc(MAXXGRID*sizeof(double)); xclon = (double *)malloc(MAXXGRID*sizeof(double)); xclat = (double *)malloc(MAXXGRID*sizeof(double)); for(i=0; i