#include #include #include #include "gridnav.h" #define MISSING -999 #define PI 3.1415927 #define RAD_TO_DEG 57.29577951 /* #define EARTH_RAD 6371.200 */ #define EARTH_RAD 6370.000 int fill_proj_parms(GridNav *gridnav); /***************************************************************************** * Fills up the gridnav structure using arguments input to the function * * input: * central_lat - the central latitude of the projection, in degrees * central_lon - the central longitude of the projection, in degrees * projection - the projection number (defined by #define's in gridnav.h) * truelat1 - the first "true" latitude. Only has an effect with * polar stereographic and lambert projections * truelat2 - the second true latitude. Only has an effect with lambert * projection * num_columns - number of columns in grid * num_rows - number of rows in grid * dx,dy - east-west (dx) and north-south (dy) distance between grid * points, in degrees for latlon (i.e., cylindrical * equidistant) projection, in km for all other projections * (including mercator). * lat_origin - latitude of grid point defined in origin_column, origin_row * lon_origin - longitude of grid point defined in origin_column, origin_row * origin_column - column for lat_origin, long_origin pair * origin_row - row for lat_origin, long_origin pair * * output: * gridnav - filled gridnav structure * *****************************************************************************/ int GRID_init(float central_lat, float central_lon, int projection, float truelat1, float truelat2, int num_columns, int num_rows, float dx, float dy, float lat_origin, float lon_origin, float origin_column, float origin_row, GridNav *gridnav) { int status = 1; gridnav->proj.central_lon = central_lon; gridnav->proj.central_lat = central_lat; gridnav->proj.map_proj = projection; gridnav->proj.truelat1 = truelat1; gridnav->proj.truelat2 = truelat2; gridnav->grid.num_columns = num_columns; gridnav->grid.num_rows = num_rows; gridnav->grid.dx = dx; gridnav->grid.dy = dy; gridnav->grid.lat_origin = lat_origin; gridnav->grid.lon_origin = lon_origin; gridnav->grid.origin_column = origin_column; gridnav->grid.origin_row = origin_row; fill_proj_parms(gridnav); return status; } /***************************************************************************** * Outputs the latitude and longitude of a grid point. * * input: * *gridnav - a pointer to a filled gridnav structure. The gridnav * structure should have been filled using one of the init * functions. * column - the column in the grid to retrieve. * row - the row in the grid to retrieve. * * output: * *lat - pointer to output latitude * *lon - pointer to output longitude * *****************************************************************************/ int GRID_to_latlon(GridNav *gridnav, float column, float row, float *lat, float *lon) { /* * Calculate the latitude and longitude for an input grid point location */ double X, Y, R; int status = 1; double eps = 0.00001; switch (gridnav->proj.map_proj) { case GRID_MERCATOR: *lat = 2 * RAD_TO_DEG * atan(exp((gridnav->proj_transform.parm2 + gridnav->grid.dx * (row - gridnav->grid.origin_row)) / gridnav->proj_transform.parm1 )) - 90; *lon = gridnav->grid.lon_origin + RAD_TO_DEG * gridnav->grid.dx * (column - gridnav->grid.origin_column) / gridnav->proj_transform.parm1; break; case GRID_LAMCON: X = (column - gridnav->grid.origin_column) * gridnav->grid.dx + gridnav->proj_transform.parm6; Y = (row - gridnav->grid.origin_row) * gridnav->grid.dy + gridnav->proj_transform.parm3 + gridnav->proj_transform.parm7; R = sqrt(X*X + Y*Y); *lat = gridnav->proj_transform.parm5 * 90 - 2 * RAD_TO_DEG * atan(gridnav->proj_transform.parm4 * pow(R, 1 / gridnav->proj_transform.parm2)); *lon = gridnav->proj.central_lon + RAD_TO_DEG / gridnav->proj_transform.parm2 * atan(X / (gridnav->proj_transform.parm5 * -Y)); while (*lon > 180) *lon -= 360; while (*lon <= -180) *lon += 360; break; case GRID_POLSTR: X = (column - gridnav->grid.origin_column) * gridnav->proj_transform.parm3 + gridnav->proj_transform.parm1; Y = (row - gridnav->grid.origin_row) * gridnav->proj_transform.parm3 + gridnav->proj_transform.parm2 + eps; *lon = gridnav->proj_transform.parm5 * -1 * atan(X / Y) * RAD_TO_DEG + gridnav->proj.central_lon; *lat = (gridnav->proj_transform.parm5 * 90) - 2 * RAD_TO_DEG * atan(X / (2 * EARTH_RAD * sin((*lon - gridnav->proj.central_lon ) / RAD_TO_DEG) + eps) ); while (*lon > 180) *lon -= 360; while (*lon <= -180) *lon += 360; break; case GRID_LATLON: *lat = (row - gridnav->grid.origin_row)*gridnav->grid.dy + gridnav->grid.lat_origin; *lon = (column - gridnav->grid.origin_column)*gridnav->grid.dx + gridnav->grid.lon_origin; break; default: /* * Unsupported map projection: set lat-lon to no-data values. */ fprintf(stderr,"GRID_to_latlon: Unsupport map projection type %d\n", gridnav->proj.map_proj); *lon = -9999; *lat = -9999; status = 0; break; } /* end of switch on map projection type */ return status; } /***************************************************************************** * Outputs grid point indices, given lat/lon location. * * input: * *gridnav - a pointer to a filled gridnav structure. The gridnav * structure should have been filled using one of the init * functions. * lat - latitude for location * lon - longitude for location * output: * *column - pointer to value for column * *row - pointer to value for row * *****************************************************************************/ int GRID_from_latlon(GridNav *gridnav, float lat, float lon, float *column, float *row) { double X, Y, Rs; double lat_rad; int status = 1; switch (gridnav->proj.map_proj) { case GRID_MERCATOR: X = gridnav->proj_transform.parm1 * ((lon - gridnav->grid.lon_origin) / RAD_TO_DEG); *column = gridnav->grid.origin_column + X / gridnav->grid.dx; lat_rad = lat / RAD_TO_DEG; Y = gridnav->proj_transform.parm1 * log((1 + sin(lat_rad)) / cos(lat_rad)); *row = gridnav->grid.origin_row + ((Y-gridnav->proj_transform.parm2) / gridnav->grid.dy);; break; case GRID_LAMCON: Rs = (EARTH_RAD / gridnav->proj_transform.parm2) * sin(gridnav->proj_transform.parm1) * pow(tan((gridnav->proj_transform.parm5 * 90 - lat) / (2 * RAD_TO_DEG))/ tan(gridnav->proj_transform.parm1 / 2), gridnav->proj_transform.parm2); *row = gridnav->grid.origin_row - (1 / gridnav->grid.dy) * (gridnav->proj_transform.parm3 + Rs * cos(gridnav->proj_transform.parm2* (lon - gridnav->proj.central_lon) / RAD_TO_DEG)) - gridnav->proj_transform.parm7 / gridnav->grid.dy; *column = gridnav->grid.origin_column + ( gridnav->proj_transform.parm5* ((Rs / gridnav->grid.dx)* sin(gridnav->proj_transform.parm2 * (lon - gridnav->proj.central_lon) / RAD_TO_DEG) - gridnav->proj_transform.parm6 / gridnav->grid.dx)); break; case GRID_POLSTR: Rs = 2 * EARTH_RAD * tan ( (45 - fabs(lat)/2) / RAD_TO_DEG ); Y = gridnav->proj_transform.parm5 * -1 * Rs * cos ((lon - gridnav->proj.central_lon) / RAD_TO_DEG); X = Rs * sin ((lon - gridnav->proj.central_lon) / RAD_TO_DEG); *row = (Y - gridnav->proj_transform.parm2) / gridnav->proj_transform.parm3 + gridnav->grid.origin_row; *column = (X - gridnav->proj_transform.parm1) / gridnav->proj_transform.parm3 + gridnav->grid.origin_column; break; case GRID_LATLON: *row = ((lat - gridnav->grid.lat_origin) / gridnav->grid.dy) + gridnav->grid.origin_row; /* If lon is negative, make it positive */ while (lon < 0) lon += 360.; *column = ((lon - gridnav->grid.lon_origin) / gridnav->grid.dx) + gridnav->grid.origin_column; break; default: fprintf(stderr,"GRID_from_latlon: Unsupported map projection type %d\n", gridnav->proj.map_proj); *column = -9999; *row = -9999; status = 0; break; } /* End of switch on map projection type */ return status; } int fill_proj_parms(GridNav *gridnav) { double orig_lat_rad; double R_orig; double r_not; double r_truelat1; int hemifactor; switch (gridnav->proj.map_proj) { case GRID_MERCATOR: gridnav->proj_transform.parm1 = EARTH_RAD * cos(gridnav->proj.truelat1 / RAD_TO_DEG); orig_lat_rad = (gridnav->grid.lat_origin) / RAD_TO_DEG; gridnav->proj_transform.parm2 = gridnav->proj_transform.parm1 * log((1. + sin(orig_lat_rad)) / cos(orig_lat_rad)); gridnav->proj_transform.parm3 = MISSING; gridnav->proj_transform.parm4 = MISSING; gridnav->proj_transform.parm5 = MISSING; break; case GRID_LAMCON: if (gridnav->proj.truelat1 >= 0) { hemifactor = 1; } else { hemifactor = -1; } /* This is Psi1 in MM5 speak */ gridnav->proj_transform.parm1 = hemifactor*(PI/2 - fabs(gridnav->proj.truelat1) / RAD_TO_DEG); /* This is Kappa in MM5 speak */ if (fabs(gridnav->proj.truelat1 - gridnav->proj.truelat2) < .01) { gridnav->proj_transform.parm2 = fabs(sin(gridnav->proj.truelat1 / RAD_TO_DEG)); } else { gridnav->proj_transform.parm2 = (log10(cos(gridnav->proj.truelat1 / RAD_TO_DEG)) - log10(cos(gridnav->proj.truelat2 / RAD_TO_DEG))) / (log10(tan((45 - fabs(gridnav->proj.truelat1) / 2) / RAD_TO_DEG)) - log10(tan((45 - fabs(gridnav->proj.truelat2) / 2) / RAD_TO_DEG))); } /* This is Yc in MM5 speak */ gridnav->proj_transform.parm3 = -EARTH_RAD/gridnav->proj_transform.parm2 * sin(gridnav->proj_transform.parm1) * pow( (tan((hemifactor * 90 - gridnav->grid.lat_origin) / (RAD_TO_DEG * 2)) / tan(gridnav->proj_transform.parm1 / 2)), gridnav->proj_transform.parm2); gridnav->proj_transform.parm4 = tan(gridnav->proj_transform.parm1 / 2)* pow(hemifactor * (gridnav->proj_transform.parm2)/ (EARTH_RAD * sin(gridnav->proj_transform.parm1)), 1 / gridnav->proj_transform.parm2); gridnav->proj_transform.parm5 = hemifactor; R_orig = (EARTH_RAD / gridnav->proj_transform.parm2) * sin(gridnav->proj_transform.parm1) * pow(tan((gridnav->proj_transform.parm5 * 90 - gridnav->grid.lat_origin) / (2 * RAD_TO_DEG))/ tan(gridnav->proj_transform.parm1 / 2), gridnav->proj_transform.parm2); /* X origin */ gridnav->proj_transform.parm6 = R_orig * sin(gridnav->proj_transform.parm2 * (gridnav->grid.lon_origin - gridnav->proj.central_lon) / RAD_TO_DEG); /* Y origin */ gridnav->proj_transform.parm7 = -(gridnav->proj_transform.parm3 + R_orig * cos(gridnav->proj_transform.parm2 * (gridnav->grid.lon_origin - gridnav->proj.central_lon) / RAD_TO_DEG)); break; case GRID_POLSTR: if (gridnav->proj.central_lat > 0) { hemifactor = 1; } else { hemifactor = -1; } /* Calculate X for origin */ r_not = 2 * EARTH_RAD * tan((45 - fabs(gridnav->grid.lat_origin) / 2) / RAD_TO_DEG); gridnav->proj_transform.parm1 = r_not * sin ( (gridnav->grid.lon_origin - gridnav->proj.central_lon) / RAD_TO_DEG); /* Calculate Y for origin */ gridnav->proj_transform.parm2 = hemifactor * -1 * r_not * cos ( (gridnav->grid.lon_origin - gridnav->proj.central_lon) / RAD_TO_DEG); /* Calculate grid spacing at pole */ r_truelat1 = 2 * EARTH_RAD * tan((45 - fabs(gridnav->proj.truelat1) / 2) / RAD_TO_DEG); gridnav->proj_transform.parm3 = gridnav->grid.dx * r_truelat1 / ( EARTH_RAD * cos (gridnav->proj.truelat1 / RAD_TO_DEG)); gridnav->proj_transform.parm4 = MISSING; gridnav->proj_transform.parm5 = hemifactor; break; case GRID_LATLON: gridnav->proj_transform.parm1 = MISSING; gridnav->proj_transform.parm2 = MISSING; gridnav->proj_transform.parm3 = MISSING; gridnav->proj_transform.parm4 = MISSING; gridnav->proj_transform.parm5 = MISSING; break; default: fprintf(stderr,"GRID_init_mm5data: Invalid Projection\n"); return 0; } return 1; } /***************************************************************************** * Rotates u and v components of wind, from being relative to earth, to * relative to grid. * * input: * *gridnav - a pointer to a filled gridnav structure. The gridnav * structure should have been filled using one of the init * functions. * lon - longitude for location * u_earth - the u component of the wind in earth (north-relative) * coordinates. * v_earth - the v component of the wind in earth (north-relative) * coordinates. * output: * *u_grid - pointer to value for u in grid coordinates * *v_grid - pointer to value for v in grid coordinates * *****************************************************************************/ int GRID_rotate_from_earth_coords(GridNav *gridnav, float lon, float u_earth, float v_earth, float *u_grid, float *v_grid) { float speed, dir; float dir_grid; /* Calculate Speed and Direction from u,v */ switch (gridnav->proj.map_proj) { case GRID_MERCATOR: *u_grid = u_earth; *v_grid = v_earth; break; case GRID_POLSTR: case GRID_LAMCON: speed = sqrt(u_earth * u_earth + v_earth * v_earth); dir = RAD_TO_DEG * atan2(-u_earth,-v_earth); while (dir >= 360) dir -= 360; while (dir < 0) dir += 360; dir_grid = dir - (lon - gridnav->proj.central_lon) * gridnav->proj_transform.parm2; while (dir_grid >= 360) dir_grid -= 360; while (dir_grid < 0) dir_grid += 360; *u_grid = -1. * speed * sin(dir_grid / RAD_TO_DEG); *v_grid = -1. * speed * cos(dir_grid / RAD_TO_DEG); break; default: fprintf(stderr, "GRID_rotate_from_earth_coords: Invalid Projection\n"); return 0; } /* End of switch projection */ return 1; } /***************************************************************************** * Rotates u and v components of wind, from being relative to earth, to * relative to grid. * * input: * *gridnav - a pointer to a filled gridnav structure. The gridnav * structure should have been filled using one of the init * functions. * lon - longitude for location * u_grid - the u component of the wind in grid coordinates * v_grid - the v component of the wind in grid coordinates * * output: * *u_earth - pointer to value for u in earth-relative coordinates * *v_grid - pointer to value for v in earth-relative coordinates * *****************************************************************************/ int GRID_rotate_to_earth_coords(GridNav *gridnav, float lon, float u_grid, float v_grid, float *u_earth, float *v_earth) { float speed, dir_grid; float dir_earth; /* Calculate Speed and Direction from u,v */ switch (gridnav->proj.map_proj) { case GRID_MERCATOR: *u_earth = u_grid; *v_earth = v_grid; break; case GRID_POLSTR: case GRID_LAMCON: speed = sqrt(u_grid * u_grid + v_grid * v_grid); dir_grid = RAD_TO_DEG * atan2(-u_grid, -v_grid); while (dir_grid >= 360) dir_grid -= 360; while (dir_grid < 0) dir_grid += 360; dir_earth = dir_grid + (lon - gridnav->proj.central_lon) * gridnav->proj_transform.parm2; while (dir_earth >= 360) dir_earth -= 360; while (dir_earth < 0) dir_earth += 360; *u_earth = -1. * speed * sin(dir_earth / RAD_TO_DEG); *v_earth = -1. * speed * cos(dir_earth / RAD_TO_DEG); break; default: fprintf(stderr, "GRID_rotate_to_earth_coords: Invalid Projection\n"); return 0; } /* End of switch projection */ return 1; }