diff --git a/regression_scripts/run_regression.sh b/regression_scripts/run_regression.sh index c9339b58..3c311357 100755 --- a/regression_scripts/run_regression.sh +++ b/regression_scripts/run_regression.sh @@ -13,7 +13,7 @@ regressiontype=$1 condaenv=$2 -if [ "$1" == "-h" ] || [ "$#" -ne 2 ] ; then +if [[ $1 == -h ]] || [[ $# != 2 ]] ; then echo "" echo "Usage: `basename $0` [-h] REGRESSIONTYPE CONDAENV" echo " Runs quick (fast, does less) or " @@ -25,16 +25,16 @@ if [ "$1" == "-h" ] || [ "$#" -ne 2 ] ; then echo " REGRESSIONTYPE: quick or daily" echo " CONDENV: name of conda env with python2 and nose" echo "" - exit 1 + return 1 fi -if [ "$regressiontype" != "quick" ] && [ "$regressiontype" != "daily" ] ; then +if [[ $regressiontype != quick ]] && [[ $regressiontype != daily ]] ; then echo "`basename $0`: invalid regressiontype, should be quick or daily" - exit 1 + return 1 fi maketarget=${regressiontype} -if [ "$regressiontype" == "daily" ] ; then +if [[ $regressiontype == daily ]] ; then maketarget="ease" fi @@ -53,7 +53,7 @@ pwd make csu_${maketarget} if [ $? -ne 0 ]; then echo "`basename $0`: csu_${maketarget} make failed" - exit -1 + return -1 fi echo "`basename $0`: running setup" cd ${PMESDR_TOP_DIR}/src/prod/meas_meta_setup @@ -61,7 +61,7 @@ pwd make csu_${maketarget} if [ $? -ne 0 ]; then echo "`basename $0`: csu_$maketarget setup failed" - exit -1 + return -1 fi echo "`basename $0`: running sir" cd ${PMESDR_TOP_DIR}/src/prod/meas_meta_sir @@ -69,10 +69,10 @@ pwd make csu_${maketarget} if [ $? -ne 0 ]; then echo "`basename $0`: csu_${maketarget} SIR failed" - exit -1 + return -1 fi make csu_${maketarget}_validate if [ $? -ne 0 ]; then echo "`basename $0`: csu_$maketarget SIR validate failed" - exit -1 + return -1 fi diff --git a/src/prod/Makefile b/src/prod/Makefile index 4942b104..5ad60172 100644 --- a/src/prod/Makefile +++ b/src/prod/Makefile @@ -66,7 +66,7 @@ CONDAENV = ${PMESDR_CONDAENV} #------------------------------------------------------------------------ SUBDIRS = calcalcs/src utils/src cetb_file/src gsx/src \ - ../sirlib/c meas_meta_make meas_meta_setup meas_meta_sir + meas_meta_make meas_meta_setup meas_meta_sir all install clean: for i in $(SUBDIRS); do \ diff --git a/src/prod/meas_meta_make/Makefile b/src/prod/meas_meta_make/Makefile index 4cad00d1..087d2133 100644 --- a/src/prod/meas_meta_make/Makefile +++ b/src/prod/meas_meta_make/Makefile @@ -13,6 +13,7 @@ # TOPDIR = $(PMESDR_TOP_DIR) BINDIR = $(TOPDIR)/bin +LIBDIR = $(TOPDIR)/lib INCDIR = $(TOPDIR)/include TESTINPUT = $(PMESDR_TOP_DIR)/src/test_inputs @@ -47,20 +48,18 @@ SYSLIBS = -lm # # end configuration section #------------------------------------------------------------------------ -MEAS_INCDIR = ../../sirlib/c/include -MEAS_LIBDIR = ../../sirlib/c -LIBS = -L$(MEAS_LIBDIR) -lcsir $(SYSLIBS) -CFLAGS = -I$(INCDIR) -I$(MEAS_INCDIR) -m64 $(CONFIG_CFLAGS) -Wall -Wno-unused-variable +LIBS = -L$(LIBDIR) -lutils $(SYSLIBS) +CFLAGS = -I$(INCDIR) -m64 $(CONFIG_CFLAGS) -Wall -Wno-unused-variable ifeq ($(LOCALE),ALPINEicc) $(info Build location: $(LOCALE)) - CFLAGS = -D'ALPINEicc=1' -I$(INCDIR) -I$(MEAS_INCDIR) $(CONFIG_CFLAGS) -O3 -ipo -xHost -Wall -Wremarks -wd981 -no-prec-div -wd2415 + CFLAGS = -D'ALPINEicc=1' -I$(INCDIR) $(CONFIG_CFLAGS) -O3 -ipo -xHost -Wall -Wremarks -wd981 -no-prec-div -wd2415 endif -DEPEND_LIBS = $(MEAS_LIBDIR)/libcsir.a +DEPEND_LIBS = $(LIBDIR)/libutils.a SRCS = meas_meta_make.c OBJS = meas_meta_make.o -HDRS = $(MEAS_INCDIR)/sir_geom.h $(INCDIR)/cetb.h +HDRS = $(INCDIR)/cetb.h $(INCDIR)/utils.h BIN = meas_meta_make all : $(BIN) @@ -68,6 +67,9 @@ all : $(BIN) $(BIN) : $(OBJS) $(HDRS) $(DEPEND_LIBS) $(CC) $(CFLAGS) -o $@ $(OBJS) $(LIBS) +meas_meta_make.o : $(SRCS) $(HDRS) + $(CC) -c $(CFLAGS) -o $@ $(SRCS) + install : $(MKDIR) $(BINDIR) $(INSTALL) $(BIN) $(BINDIR) diff --git a/src/prod/meas_meta_make/meas_meta_make.c b/src/prod/meas_meta_make/meas_meta_make.c index b07dd8f4..146bcade 100644 --- a/src/prod/meas_meta_make/meas_meta_make.c +++ b/src/prod/meas_meta_make/meas_meta_make.c @@ -11,15 +11,11 @@ #include #include #include -#ifdef JANUSicc -#include -#else #include -#endif #include -#include -#include +#include "cetb.h" +#include "utils.h" #define prog_version 1.2 /* program version */ #define prog_name "meas_meta_make" @@ -32,13 +28,6 @@ #define TRUE 1 #define FALSE 0 -/****************************************************************************/ -/* default location of the SIR standard region definition */ - -char rname[] = "regiondef1.dat"; /* file defining region codes */ - -/********************************************************************/ - /* function prototypes */ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, @@ -46,7 +35,7 @@ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, static int get_file_names(FILE *mout, int *argn, char *argv[]); -static void get_region_data(int regnum, int resolution_ind, int *iproj, int *dateline, float *latl, +static void get_region_data(int regnum, int resolution_ind, int *iproj, float *latl, float *lonl, float *lath, float *lonh, char *regname); /****************************************************************************/ @@ -254,7 +243,6 @@ static int get_file_names(FILE *mout, int *argn, char *argv[]) * * output: * int *iproj - short version of region id (i.e. 8 rather than 308 etc) - * int dateline - this is always zero for our projections - might not be needed * float *latl, *lonl - latitude and longitude of lower left corner * of lower left corner pixel * float *lath, *lonh - latitude and longitude of upper right corner @@ -262,7 +250,7 @@ static int get_file_names(FILE *mout, int *argn, char *argv[]) * char *regname - EASE2_N etc. */ static void get_region_data(int regnum, int resolution_ind, int *iproj, - int *dateline, float *latl, float *lonl, + float *latl, float *lonl, float *lath, float *lonh, char *regname) { cetb_region_id cetb_region; @@ -271,14 +259,13 @@ static void get_region_data(int regnum, int resolution_ind, int *iproj, projection_offset = regnum - CETB_PROJECTION_BASE_NUMBER; cetb_region = (cetb_region_id)( projection_offset + (CETB_NUMBER_BASE_RESOLUTIONS * resolution_ind) ); - dateline = 0; strncpy( regname, cetb_region_id_name[ resolution_ind ][ projection_offset ], CETB_MAX_LEN_REGION_ID_NAME ); - *latl = cetb_latitude_extent[ resolution_ind ][ projection_offset ][0]; - *lath = cetb_latitude_extent[ resolution_ind ][ projection_offset ][1]; - *lonl = cetb_longitude_extent[ resolution_ind ][ projection_offset ][0]; - *lonh = cetb_longitude_extent[ resolution_ind ][ projection_offset ][1]; + *latl = (float)cetb_latitude_extent[ resolution_ind ][ projection_offset ][0]; + *lath = (float)cetb_latitude_extent[ resolution_ind ][ projection_offset ][1]; + *lonl = (float)cetb_longitude_extent[ resolution_ind ][ projection_offset ][0]; + *lonh = (float)cetb_longitude_extent[ resolution_ind ][ projection_offset ][1]; switch ( cetb_region ) { case CETB_EASE2_N: @@ -338,11 +325,10 @@ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, char rfile[250]; int pfile; FILE *pid; - int nregions, poleflag, dateline, iproj, regnum, iregion, ircnt; + int nregions, poleflag, iproj, regnum, iregion, ircnt; float latl, lonl, lath, lonh; int projt, nsx, nsy, xdim, ydim, nease; float ascale, bscale, xdeg, ydeg, a0, b0, aorglon, aorglat; - float maplxlon, maprxlon; int iasc, ibeam, ipolar; int non_size; int nsect, nsx2, nsy2; @@ -426,7 +412,6 @@ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, ircnt=0; /* count total regions */ for (iregion=0; iregion 0) { /* get region definition from contents of cetb.h */ - get_region_data( regnum, resolution_ind, &iproj, &dateline, &latl, &lonl, + get_region_data( regnum, resolution_ind, &iproj, &latl, &lonl, &lath, &lonh, regname ); if (((regnum >= 0) && (regnum < 100)) || (regnum >= 120)) poleflag=0; /* non-polar area */ - fprintf( stderr, "%s: Region name: '%s' Def Proj %d Dateline %d\n", - __FUNCTION__, regname, iproj, dateline ); + fprintf( stderr, "%s: Region name: '%s' Def Proj %d \n", + __FUNCTION__, regname, iproj ); } else { fprintf( stderr, "%s: Region is not defined and out of bounds\n", __FUNCTION__ ); exit(-1); @@ -452,15 +437,6 @@ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, fprintf( stderr, "%s: Longitude range: %f %f\n", __FUNCTION__, lonl, lonh ); fprintf( stderr, "%s: Region polar code (0=arbitrary, 1=N pol, 2=S pole: %d\n", __FUNCTION__, poleflag ); - if (dateline) { - fprintf( stderr, "%s: Region crosses dateline\n", __FUNCTION__ ); - maplxlon=min(lonh,lonl); - maprxlon=max(lonh,lonl); - lonl=maplxlon; - lonh=maprxlon; - fprintf( stderr, "%s: Corrected longitude range: %f %f\n", - __FUNCTION__, lonl, lonh ); - } /* write region ID number and bound box info to meta file */ fprintf(mout," Begin_region_description\n"); @@ -469,7 +445,6 @@ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, fprintf(mout," Latitude_high=%16.9f\n", lath); fprintf(mout," Longitude_low=%16.9f\n", lonl); fprintf(mout," Longitude_high=%16.9f\n", lonh); - fprintf(mout," Dateline_crossing=%c\n", TF[dateline]); fprintf(mout," Polar_flag=%c\n", TF[poleflag]); fprintf(mout," Region_name=%s\n", regname); @@ -499,7 +474,7 @@ static int get_region_parms( FILE *mout, int *argn, char *argv[], int F_num, ind=resolution_ind; /* standard base resolution */ fprintf( stderr, "%s: EASE2 parameters: proj=%d nease=%d ind=%d\n", __FUNCTION__, projt, nease, ind); - ease2_map_info(projt, nease, ind, &map_equatorial_radius_m, + utils_ease2_map_info(projt, nease, ind, &map_equatorial_radius_m, &map_eccentricity, &e2, &map_reference_latitude, &map_reference_longitude, &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, diff --git a/src/prod/meas_meta_setup/Makefile b/src/prod/meas_meta_setup/Makefile index 150301bd..c5c28988 100644 --- a/src/prod/meas_meta_setup/Makefile +++ b/src/prod/meas_meta_setup/Makefile @@ -81,17 +81,15 @@ endif # # end configuration section #------------------------------------------------------------------------ -MEAS_INCDIR = ../../sirlib/c/include -MEAS_LIBDIR = LIBS = -L$(LIBDIR) -L$(NETCDF4_LIBDIR) -L$(HDF5_LIBDIR) -L$(UDUNITS_LIBDIR) \ - -lgsx -lutils -lcalcalcs -ludunits2 -lcsir \ + -lgsx -lutils -lcalcalcs -ludunits2 \ -lnetcdf -l$(HDFLIB1) -l$(HDFLIB2) -lz $(SYSLIBS) CFLAGS = -I$(INCDIR) -I$(NETCDF4_INCDIR) -I$(HDF5_INCDIR) -I$(UDUNITS_INCDIR) \ - -I$(MEAS_INCDIR) $(CONFIG_CFLAGS) -Wall -Wno-unused-variable + $(CONFIG_CFLAGS) -Wall -Wno-unused-variable ifeq ($(LOCALE),ALPINEicc) $(info Build location: $(LOCALE)) CFLAGS = -D'ALPINEicc=1' -I$(INCDIR) -I$(NETCDF4_INCDIR) -I$(HDF5_INCDIR) \ - -I$(UDUNITS_INCDIR) -I$(MEAS_INCDIR) -ipo -xHost -no-prec-div \ + -I$(UDUNITS_INCDIR) -ipo -xHost -no-prec-div \ -Wall -Wremarks -wd981 -wd2415 -inline-level=0 endif @@ -99,11 +97,11 @@ endif DEPEND_LIBS = $(LIBDIR)/libutils.a \ $(LIBDIR)/libcalcalcs.a \ $(UDUNITS_LIBDIR)/libudunits2.a \ - $(LIBDIR)/libgsx.a $(LIBDIR)/libcsir.a + $(LIBDIR)/libgsx.a SRCS = meas_meta_setup.c OBJS = meas_meta_setup.o -HDRS = $(INCDIR)/gsx.h $(INCDIR)/cetb.h $(MEAS_INCDIR)/sir_geom.h \ +HDRS = $(INCDIR)/gsx.h $(INCDIR)/cetb.h $(INCDIR)/utils.h \ $(INCDIR)/utCalendar2_cal.h BIN = meas_meta_setup diff --git a/src/prod/meas_meta_setup/meas_meta_setup.c b/src/prod/meas_meta_setup/meas_meta_setup.c index 0109a0b3..db891ee7 100644 --- a/src/prod/meas_meta_setup/meas_meta_setup.c +++ b/src/prod/meas_meta_setup/meas_meta_setup.c @@ -27,7 +27,6 @@ #include "utils.h" #include "gsx.h" #include "utils.h" -#include "sir_geom.h" #include "utCalendar2_cal.h" #define prog_version 0.3 /* program version */ @@ -41,8 +40,6 @@ #define RESPONSEMULT 1000 /* response multiplier */ #define HASAZIMUTHANGLE 1 /* include azimuth angle in output setup file if 1, set to 0 to not include az ang (smaller file) */ -#define DTR ((2.0*(M_PI))/360.0) /* degrees to radians */ - #define AEARTH 6378.1363 /* SEMI-MAJOR AXIS OF EARTH, a, KM */ #define FLAT 3.3528131778969144e-3 /* = 1/298.257 FLATNESS, f, f=1-sqrt(1-e**2) */ @@ -82,6 +79,14 @@ static void no_trailing_blanks(char *s) return; } +static int intfix(float r) +{ + int ret_val = r; + if (ret_val - r > 0.5) ret_val--; + if (r - ret_val > 0.5) ret_val++; + return(ret_val); +} + /****************************************************************************/ typedef enum { UNKNOWN_LTOD=-1, @@ -166,6 +171,16 @@ static int day_offset_from( int year, int month, int day, static int ltod_day_offset( int dstart, int dend, int *midDay, int *startDayOffset, int *endDayOffset, int *imageDayOffset ); +static void latlon2pix(float alon, float alat, float *x, float *y, + int iopt, float ascale, float bscale, float a0, float b0); +static void pixtolatlon(float x, float y, float *alon, float *alat, + int iopt, float ascale, float bscale, float a0, float b0); +static void ease2grid(int iopt, float alon, float alat, + float *thelon, float *thelat, float ascale, float bscale); +static void iease2grid(int iopt, float *alon, float *alat, + float thelon, float thelat, float ascale, float bscale); +static double easeconv_normalize_degrees(double b); +static void f2ipix(float x, float y, int *ix, int *iy, int nsx, int nsy); /****************************************************************************/ @@ -213,7 +228,7 @@ int main(int argc,char *argv[]) int ibeam,icc,icmax,icmin,nfile; float cx,cy,lath,latl,lonh,lonl; int nsx,nsy,projt; - float ascale,bscale,a0,b0,xdeg,ydeg,x,y; + float ascale,bscale,a0,b0,x,y; float tbmin=1.e10,tbmax=-1.e10; int iadd1, box_size, box_size_flag=0;; @@ -1067,12 +1082,10 @@ int main(int argc,char *argv[]) bscale=save_area.sav_bscale[iregion]; a0=save_area.sav_a0[iregion]; b0=save_area.sav_b0[iregion]; - ydeg=save_area.sav_ydeg[iregion]; - xdeg=save_area.sav_xdeg[iregion]; dscale=save_area.sav_km[iregion]; /* transform center lat/lon location of measurement to image pixel location */ - latlon2pix(cx, cy, &x, &y, projt, xdeg, ydeg, ascale, bscale, a0, b0); + latlon2pix(cx, cy, &x, &y, projt, ascale, bscale, a0, b0); /* quantize pixel location to 1-based integer pixel indices */ f2ipix(x, y, &ix2, &iy2, nsx, nsy); /* note: ix2,iy2 are 1-based pixel address of center */ @@ -1089,7 +1102,7 @@ int main(int argc,char *argv[]) centers to the center of the output pixel */ x=(ix2+0.5f); y=(iy2+0.5f); - pixtolatlon(x, y, &clon, &clat, projt, xdeg, ydeg, ascale, bscale, a0, b0); + pixtolatlon(x, y, &clon, &clat, projt, ascale, bscale, a0, b0); /* define size of box centered at(ix2,iy2) in which the gain response is computed for each pixel in the box and tested to see if @@ -2047,7 +2060,6 @@ void compute_locations(region_save *a, int *nregions, int **noffset, short int * for (ix=0; ixsav_nsx[iregion]; ix++) { x=(ix+1.5f); /* center of pixel, 1-based */ pixtolatlon(x, y, &clon, &clat, a->sav_projt[iregion], - a->sav_xdeg[iregion], a->sav_ydeg[iregion], a->sav_ascale[iregion], a->sav_bscale[iregion], a->sav_a0[iregion], a->sav_b0[iregion]); iadd=a->sav_nsx[iregion]*iy+ix; /* zero-based lexicographic pixel address */ @@ -2100,8 +2112,8 @@ float gsx_antenna_response(float x_rel, float y_rel, float theta, float semimajo float x, y, cross_beam_size, along_beam_size, t1, t2, weight; /* rotate coordinate system to align with look direction */ - x=(float) ( ( (cos(theta*DTR) ) * x_rel ) - ( (sin(theta*DTR) ) * y_rel ) ); - y=(float) ( ( (sin(theta*DTR) ) * x_rel ) + ( (cos(theta*DTR) ) * y_rel ) ); + x=(float) ( ( (cos(theta*UTILS_DTR) ) * x_rel ) - ( (sin(theta*UTILS_DTR) ) * y_rel ) ); + y=(float) ( ( (sin(theta*UTILS_DTR) ) * x_rel ) + ( (cos(theta*UTILS_DTR) ) * y_rel ) ); /* compute approximate antenna response Antenna weighting is estimation from SSMI Users Guide 21-27 */ @@ -2131,7 +2143,7 @@ void rel_latlon(float *x_rel, float *y_rel, float alon, float alat, float rlon, float r,r2,rel_rlat,rel_rlon; - r=(float) ( ( 1.0 - ( ( (sin(rlat*DTR) ) * (sin(rlat*DTR) ) ) * FLAT ) ) * AEARTH ); + r=(float) ( ( 1.0 - ( ( (sin(rlat*UTILS_DTR) ) * (sin(rlat*UTILS_DTR) ) ) * FLAT ) ) * AEARTH ); rel_rlat=alat-rlat; rel_rlon=alon-rlon; @@ -2141,9 +2153,9 @@ void rel_latlon(float *x_rel, float *y_rel, float alon, float alat, float rlon, else rel_rlon=rel_rlon+360.f; } - r2=r*(float)(cos(rlat*DTR)); - *x_rel=(float)(r2*(sin(rel_rlon*DTR))); - *y_rel=(float)((r*(sin(rel_rlat*DTR)))+((1.-(cos(rel_rlon*DTR)))*(sin(rlat*DTR))*r2)); + r2=r*(float)(cos(rlat*UTILS_DTR)); + *x_rel=(float)(r2*(sin(rel_rlon*UTILS_DTR))); + *y_rel=(float)((r*(sin(rel_rlat*UTILS_DTR)))+((1.-(cos(rel_rlon*UTILS_DTR)))*(sin(rlat*UTILS_DTR))*r2)); } /* *********************************************************************** */ @@ -2188,7 +2200,7 @@ float km2pix(float *x, float *y, int iopt, float ascale, float bscale, int *stat case 10: /* EASE2 T */ nease=ascale; ind=bscale; - ease2_map_info(iopt, nease, ind, &map_equatorial_radius_m, + utils_ease2_map_info(iopt, nease, ind, &map_equatorial_radius_m, &map_eccentricity, &e2, &map_reference_latitude, &map_reference_longitude, &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, @@ -3152,7 +3164,7 @@ static int ltod_day_offset( int dstart, int dend, int *midDay, int *startDayOffset, int *endDayOffset, int *imageDayOffset ) { - *midDay = dstart + round((dend-dstart)/2); + *midDay = dstart + (int)round((dend-dstart)/2); if ( dstart == dend ) { *startDayOffset = -1; @@ -3290,3 +3302,265 @@ static int check_for_consistent_regions( region_save *save_area, } +static void pixtolatlon(float x, float y, float *alon, float *alat, + int iopt, float ascale, float bscale, float a0, float b0) +{ + /* computes the lat, long position of the lower-left corner of the + x,y-th pixel. pixels indexed [1..nsx,1..nsy] */ + + float thelon, thelat; + + switch(iopt) { + case 8: + case 9: + case 10: + thelon = x - (float)1.0 + a0; + thelat = y - (float)1.0 + b0; + iease2grid(iopt, alon, alat, thelon, thelat, ascale, bscale); + break; + default: + *alon=0.0; + *alat=0.0; + } + return; +} + +static void latlon2pix(float alon, float alat, float *x, float *y, + int iopt, float ascale, float bscale, float a0, float b0) +{ + /* computes the x,y pixel coordinates given the lat, lon position + x,y are fractional values not limited to within image */ + + static float thelon, thelat; + + switch(iopt) { + case 8: + case 9: + case 10: + ease2grid(iopt, alon, alat, &thelon, &thelat, ascale, bscale); + *x = thelon + (float)1.0 - a0; + *y = thelat + (float)1.0 - b0; + break; + default: + *x=0.0; + *y=0.0; + } + return; +} + +static void ease2grid(int iopt, float alon, float alat, + float *thelon, float *thelat, float ascale, float bscale) +{ +/* + COMPUTE THE FORWARD "EASE2" GRID TRANSFORMATION + + GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND + THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED + USING THE "EASE2 GRID" (VERSION 2.0) TRANSFORMATION GIVEN IN IDL + SOURCE CODE SUPPLIED BY MJ BRODZIK + RADIUS EARTH=6378.137 KM (WGS 84) + MAP ECCENTRICITY=0.081819190843 (WGS84) + + inputs: + iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M + alon, alat: lon, lat (deg) to convert (can be outside of image) + ascale and bscale should be integer valued) + ascale: grid scale factor (0..5) pixel size is (bscale/2^ascale) + bscale: base grid scale index (ind=int(bscale)) + + outputs: + thelon: X coordinate in pixels (can be outside of image) + thelat: Y coordinate in pixels (can be outside of image) + +*/ + double map_equatorial_radius_m,map_eccentricity, e2, + map_reference_latitude, map_reference_longitude, + map_second_reference_latitude, sin_phi1, cos_phi1, kz, + map_scale, r0, s0, epsilon; + int bcols, brows; + + int ind = intfix(bscale); + int isc = intfix(ascale); + double dlon = alon; + double phi = UTILS_DTR * alat; + double lam = dlon; + + double sin_phi, q, qp, rho, x, y; + + /* get base EASE2 map projection parameters */ + utils_ease2_map_info(iopt, isc, ind, &map_equatorial_radius_m, &map_eccentricity, + &e2, &map_reference_latitude, &map_reference_longitude, + &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, + &map_scale, &bcols, &brows, &r0, &s0, &epsilon); + + dlon = dlon - map_reference_longitude; + dlon = easeconv_normalize_degrees( dlon ); + lam = UTILS_DTR * dlon; + + sin_phi=sin(phi); + q = ( 1.0 - e2 ) * ( ( sin_phi / ( 1.0 - e2 * sin_phi * sin_phi ) ) + - ( 1.0 / ( 2.0 * map_eccentricity ) ) + * log( ( 1.0 - map_eccentricity * sin_phi ) + / ( 1.0 + map_eccentricity * sin_phi ) ) ); + + switch (iopt) { + case 8: /* EASE2 grid north */ + qp = 1.0 - ( ( 1.0 - e2 ) / ( 2.0 * map_eccentricity ) + * log( ( 1.0 - map_eccentricity ) + / ( 1.0 + map_eccentricity ) ) ); + if ( abs( qp - q ) < epsilon ) + rho = 0.0; + else + rho = map_equatorial_radius_m * sqrt( qp - q ); + x = rho * sin( lam ); + y = -rho * cos( lam ); + break; + case 9: /* EASE2 grid south */ + qp = 1.0 - ( ( 1.0 - e2 ) / ( 2.0 * map_eccentricity ) + * log( ( 1.0 - map_eccentricity ) + / ( 1.0 + map_eccentricity ) ) ); + if ( abs( qp + q ) < epsilon ) + rho = 0.0; + else + rho = map_equatorial_radius_m * sqrt( qp + q ); + x = rho * sin( lam ); + y = rho * cos( lam ); + break; + case 10: /* EASE2 cylindrical */ + x = map_equatorial_radius_m * kz * lam; + y = ( map_equatorial_radius_m * q ) / ( 2.0 * kz ); + break; + default: + fprintf(stderr,"*** invalid EASE2 projection specificaion %d in ease2grid\n",iopt); + break; + } + + *thelon = (float) (r0 + ( x / map_scale ) + 0.5); + *thelat = (float) (s0 + ( y / map_scale ) + 0.5); /* 0 at bottom (BYU SIR files) */ + + return; +} + +static void iease2grid(int iopt, float *alon, float *alat, + float thelon, float thelat, float ascale, float bscale) +{ +/* + COMPUTE THE INVERSE "EASE2" GRID TRANSFORM + + GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND + THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED + USING THE "EASE GRID" (VERSION 2.0) TRANSFORMATION GIVEN IN IDL + SOURCE CODE SUPPLIED BY MJ BRODZIK + + inputs: + iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M + thelon: X coordinate in pixels (can be outside of image) + thelat: Y coordinate in pixels (can be outside of image) + ascale and bscale should be integer valued) + ascale: grid scale factor (0..5) pixel size is (bscale/2^ascale) + bscale: base grid scale index (ind=int(bscale)) + + outputs: + alon, alat: lon, lat location in deg (can be outside of image) +*/ + double map_equatorial_radius_m,map_eccentricity, e2, + map_reference_latitude, map_reference_longitude, + map_second_reference_latitude, sin_phi1, cos_phi1, kz, + map_scale, r0, s0, epsilon; + int bcols, brows; + + int ind = intfix(bscale); + int isc = intfix(ascale); + + double lam, arg, phi, beta, qp, rho2, x, y, e4, e6; + + /* get base EASE2 map projection parameters */ + utils_ease2_map_info(iopt, isc, ind, &map_equatorial_radius_m, &map_eccentricity, + &e2, &map_reference_latitude, &map_reference_longitude, + &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, + &map_scale, &bcols, &brows, &r0, &s0, &epsilon); + e4 = e2 * e2; + e6 = e4 * e2; + + /* qp is the function q evaluated at phi = 90.0 deg */ + qp = ( 1.0 - e2 ) * ( ( 1.0 / ( 1.0 - e2 ) ) + - ( 1.0 / ( 2.0 * map_eccentricity ) ) + * log( ( 1.0 - map_eccentricity ) + / ( 1.0 + map_eccentricity ) ) ); + + x = ((double) thelon - r0 - 0.5) * map_scale; + //y = (s0 - (double) thelat + 0.5) * map_scale; /* 0 at top (NSIDC) */ + y = ((double) thelat - 0.5 - s0) * map_scale; /* 0 at bottom (BYU SIR files) */ + + switch (iopt) { + case 8: /* EASE2 grid north */ + rho2 = ( x * x ) + ( y * y ); + arg=1.0 - ( rho2 / ( map_equatorial_radius_m * map_equatorial_radius_m * qp ) ); + if (arg > 1.0) arg=1.0; + if (arg < -1.0) arg=-1.0; + beta = asin( arg ); + lam = atan2( x, -y ); + break; + case 9: /* EASE2 grid south */ + rho2 = ( x * x ) + ( y * y ); + arg = 1.0 - ( rho2 / ( map_equatorial_radius_m * map_equatorial_radius_m * qp ) ); + if (arg > 1.0) arg=1.0; + if (arg < -1.0) arg=-1.0; + beta = -asin( arg ); + lam = atan2( x, y ); + break; + case 10: /* EASE2 cylindrical */ + arg = 2.0 * y * kz / ( map_equatorial_radius_m * qp ); + if (arg > 1.0) arg=1.0; + if (arg < -1.0) arg=-1.0; + beta = asin( arg ); + lam = x / ( map_equatorial_radius_m * kz ); + break; + default: + fprintf(stderr,"*** invalid EASE2 projection specification %d in iease2grid\n",iopt); + break; + } + + phi = beta + + ( ( ( e2 / 3.0 ) + ( ( 31.0 / 180.0 ) * e4 ) + + ( ( 517.0 / 5040.0 ) * e6 ) ) * sin( 2.0 * beta ) ) + + ( ( ( ( 23.0 / 360.0 ) * e4) + + ( ( 251.0 / 3780.0 ) * e6 ) ) * sin( 4.0 * beta ) ) + + ( ( ( 761.0 / 45360.0 ) * e6 ) * sin( 6.0 * beta ) ); + + *alat = (float) (UTILS_RTD * phi); + *alon = (float) (easeconv_normalize_degrees( map_reference_longitude + ( UTILS_RTD*lam ) ) ); + + return; +} + +static double easeconv_normalize_degrees(double b) +{ /* return -180 <= b <= 180.0 */ + while (b < -180.0) + b = b + 360.0; + while (b > 180.0) + b = b - 360.0; + return(b); +} + +/* floating point to integer pixel location quantization routine */ + +static void f2ipix(float x, float y, int *ix, int *iy, int nsx, int nsy) +{ + /* quantizes the floating point pixel location to the actual pixel value + returns a zero if location is outside of image limits + a small amount (0.002 pixels) of rounding is permitted*/ + + if (x+0.0002 >= 1.0 && x+0.0002 <= (float) (nsx+1)) + *ix = (int)floor(x+0.0002); + else + *ix = 0; + + if (y+0.0002 >= 1.0 && y+0.0002 <= (float) (nsy+1)) + *iy = (int)floor(y+0.0002); + else + *iy = 0; + + return; +} + diff --git a/src/prod/meas_meta_sir/meas_meta_sir.c b/src/prod/meas_meta_sir/meas_meta_sir.c index b3046ecd..ee40a65d 100644 --- a/src/prod/meas_meta_sir/meas_meta_sir.c +++ b/src/prod/meas_meta_sir/meas_meta_sir.c @@ -11,11 +11,7 @@ #include #include #include -#ifdef JANUSicc -#include -#else #include -#endif #include #include #include @@ -147,8 +143,6 @@ int main(int argc, char **argv) cetb_sensor_id sensor_id; cetb_direction_id direction_id=CETB_NO_DIRECTION; unsigned short tb_fill_value=CETB_NCATTS_TB_FILL_VALUE; - unsigned short stokes_fill_value=CETB_NCATTS_STOKES_FILL_VALUE; - unsigned short tb_or_stokes_fill_value; unsigned short tb_valid_range[ 2 ] = { CETB_NCATTS_TB_MIN, CETB_NCATTS_TB_MAX }; unsigned short stokes_valid_range[ 2 ] = { CETB_NCATTS_STOKES_MIN, CETB_NCATTS_STOKES_MAX }; unsigned short tb_or_stokes_valid_range[ 2 ]; @@ -161,7 +155,6 @@ int main(int argc, char **argv) CETB_NCATTS_TB_NUM_SAMPLES_MAX }; short theta_fill_value=CETB_NCATTS_THETA_FILL_VALUE; short theta_valid_range[ 2 ] = { CETB_NCATTS_THETA_MIN, CETB_NCATTS_THETA_MAX }; - float error_valid_range[ 2 ] = { 0.0, NC_MAX_FLOAT }; float tb_or_stokes_scaled_min, tb_or_stokes_scaled_max, tb_or_stokes_scale_factor; float tb_or_stokes_stddev_scale_factor, tb_or_stokes_SIR_offset; int tb_or_stokes_add_offset, tb_or_stokes_stddev_add_offset; @@ -174,7 +167,7 @@ int main(int argc, char **argv) int median_flag = 0; /* default: no median filter in SIRF algorithm */ int ibeam = 0; int resolution; - cetb_resolution_id base_resolution = 0; + cetb_resolution_id base_resolution = (cetb_resolution_id)0; /* begin program */ diff --git a/src/prod/utils/src/Makefile b/src/prod/utils/src/Makefile index afb5058b..8e9f3492 100644 --- a/src/prod/utils/src/Makefile +++ b/src/prod/utils/src/Makefile @@ -69,24 +69,23 @@ endif # end configuration section #------------------------------------------------------------------------ CFLAGS = -I$(INCDIR) -m64 $(CONFIG_CFLAGS) -Wall -Wno-unused-variable +LIBS = $(SYSLIBS) ARFLAGS = ruv ifeq ($(LOCALE),ALPINEicc) $(info Build location: $(LOCALE)) CFLAGS = -D'ALPINEicc=1' -I$(INCDIR) -I$(NETCDF4_INCDIR) -I$(UDUNITS_INCDIR) $(CONFIG_CFLAGS) -O3 -ipo -g -xHost -no-prec-div -Wall -Wremarks -wd981 -wd2415 ARFLAGS = crs -endif - -DEPEND_LIBS = +endif SRCS = utils.c OBJS = utils.o HDRS = utils.h -all : libutils.a $(LIBDIR)/libutils.a +all : libutils.a $(LIBDIR)/libutils.a $(LIBDIR)/libutils.a : libutils.a -libutils.a : $(OBJS) $(DEPEND_LIBS) $(HDRS) +libutils.a : $(OBJS) $(LIBS) $(HDRS) $(AR) $(ARFLAGS) libutils.a $(OBJS) $(MKDIR) $(LIBDIR) $(INCDIR) $(INSTALL) libutils.a $(LIBDIR) diff --git a/src/prod/utils/src/utils.c b/src/prod/utils/src/utils.c index 1bf52704..60b7802d 100644 --- a/src/prod/utils/src/utils.c +++ b/src/prod/utils/src/utils.c @@ -6,6 +6,9 @@ */ #include "cetb.h" #include "utils.h" +#include + +static int dceil( double r ); /* * utils_allocate_clean_aligned_memory - allocate aligned memory that is @@ -30,3 +33,166 @@ int utils_allocate_clean_aligned_memory( void **this, size_t size ) { return 0; } + +void utils_ease2_map_info(int iopt, int isc, int ind, + double *map_equatorial_radius_m, double *map_eccentricity, + double *e2, double *map_reference_latitude, + double *map_reference_longitude, + double *map_second_reference_latitude,double * sin_phi1, + double *cos_phi1, double *kz, + double *map_scale, int *bcols, int *brows, + double *r0, double *s0, double *epsilon) +{ /* define key EASE2 grid information + + inputs + iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M + isc: scale factor 0..5 grid size is (basesize(ind))/2^isc + ind: base grid size index (map units per cell in m + + NSIDC .grd file for isc=0 + project type ind=0 ind=1 ind=2 + N EASE2_N25km EASE2_N36km EASE2_N24km + S EASE2_S25km EASE2_S36km EASE2_S24km + T/M EASE2_T25km EASE2_M36km EASE2_M24km + + cell size (m) for isc=0 (scale is reduced by 2^isc) + project type ind=0 ind=1 ind=2 + N 25000.0 36000.0 24000.0 + S 25000.0 36000.0 24000.0 + T/M T25025.26 M36032.220840584 M24021.4805603 + + for a given base cell size isc is related to NSIDC .grd file names + isc N .grd name S .grd name T .grd name + 0 EASE2_N25km EASE2_S25km EASE2_T25km + 1 EASE2_N12.5km EASE2_S12.5km EASE2_T12.5km + 2 EASE2_N6.25km EASE2_S6.25km EASE2_T6.25km + 3 EASE2_N3.125km EASE2_S3.125km EASE2_T3.125km + 4 EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km + + outputs + map_equatorial_radius_m EASE2 Earth equitorial radius (km) [WGS84] + map_eccentricity EASE2 Earth eccentricity [WGS84] + map_reference_latitude Reference latitude (deg) + map_reference_longitude Reference longitude (deg) + map_second_reference_latitude Secondary reference longitude* (deg) + sin_phi1, cos_phi1 kz EASE2 Cylin parameters* + map_scale EASE2 map projection pixel size (km) + bcols, brows, EASE2 grid size in pixels + r0, s0 EASE2 base projection size in pixels + epsilon EASE2 near-polar test factor + + *these parameters only assigned values if projection is T + + */ + + double base; + int m, nx, ny; + int region_index; + + *map_equatorial_radius_m = 6378137.0 ; /* WGS84 */ + *map_eccentricity = 0.081819190843 ; /* WGS84 */ + *e2 = *map_eccentricity * *map_eccentricity; + *map_reference_longitude = 0.0; + *epsilon = 1.e-6; + + /* map-specific parameters - these all come from cetb.h */ + switch (iopt) { + case 8: /* EASE2 grid north */ + *map_reference_latitude = 90.0; + switch(ind) { + case CETB_36KM: /* EASE2_N36km.gpd */ + region_index = CETB_36KM*CETB_NUMBER_PROJECTIONS; + base=cetb_exact_scale_m[region_index][0]; //36000.0; + nx=cetb_grid_cols[region_index][0]; //500; + ny=cetb_grid_rows[region_index][0]; //500; + break; + case CETB_24KM: /* EASE2_N24km.gpd */ + region_index = CETB_24KM*CETB_NUMBER_PROJECTIONS; + base=cetb_exact_scale_m[region_index][0]; //24000.0; + nx=cetb_grid_cols[region_index][0]; //750; + ny=cetb_grid_rows[region_index][0]; //750; + break; + default: /* EASE2_N25km.gpd */ + region_index = CETB_25KM*CETB_NUMBER_PROJECTIONS; + base=cetb_exact_scale_m[region_index][0]; //25000.0; + nx=cetb_grid_cols[region_index][0]; //720; + ny=cetb_grid_rows[region_index][0]; //720; + } + break; + case 9: /* EASE2 grid south */ + *map_reference_latitude = -90.0; + switch(ind) { + case CETB_36KM: /* EASE2_S36km.gpd */ + region_index = (CETB_36KM*CETB_NUMBER_PROJECTIONS)+1; + base=cetb_exact_scale_m[region_index][0]; //36000.0; + nx=cetb_grid_cols[region_index][0]; //500; + ny=cetb_grid_rows[region_index][0]; //500; + break; + case CETB_24KM: /* EASE2_S24km.gpd */ + region_index = (CETB_24KM*CETB_NUMBER_PROJECTIONS)+1; + base=cetb_exact_scale_m[region_index][0]; //24000.0; + nx=cetb_grid_cols[region_index][0]; //750; + ny=cetb_grid_rows[region_index][0]; //750; + break; + default: /* EASE2_S25km.gpd */ + region_index = (CETB_25KM*CETB_NUMBER_PROJECTIONS)+1; + base=cetb_exact_scale_m[region_index][0]; //25000.0; + nx=cetb_grid_cols[region_index][0]; //720; + ny=cetb_grid_rows[region_index][0]; //720; + } + break; + case 10: /* EASE2 cylindrical */ + *map_reference_latitude = 0.0; + *map_second_reference_latitude = 30.0; + *sin_phi1 = sin( UTILS_DTR * *map_second_reference_latitude ); + *cos_phi1 = cos( UTILS_DTR * *map_second_reference_latitude ); + *kz = *cos_phi1 / sqrt( 1.0 - *e2 * *sin_phi1 * *sin_phi1 ); + switch(ind) { + case CETB_36KM: /* EASE2_M36km.gpd */ + region_index = (CETB_36KM*CETB_NUMBER_PROJECTIONS)+2; + base=cetb_exact_scale_m[region_index][0]; //36032.220840584; + nx=cetb_grid_cols[region_index][0]; //964; + ny=cetb_grid_rows[region_index][0]; //406; + break; + case CETB_24KM: /* EASE2_M24km.gpd */ + region_index = (CETB_24KM*CETB_NUMBER_PROJECTIONS)+2; + base=cetb_exact_scale_m[region_index][0]; //24021.480560389347; + nx=cetb_grid_cols[region_index][0]; //1446; + ny=cetb_grid_rows[region_index][0]; //609; + break; + default: /* EASE2_T25km.gpd */ + region_index = (CETB_25KM*CETB_NUMBER_PROJECTIONS)+2; + base=cetb_exact_scale_m[region_index][0]; //25025.26000; + nx=cetb_grid_cols[region_index][0]; //1388; + ny=cetb_grid_rows[region_index][0]; //540; /* originally was 538 */ + } + } + + + /* grid info */ + if (isc>=0) { + for (m=1; isc>0; isc--) m *= 2; /* compute power-law scale factor */ + *map_scale = base / (double) m; + *bcols = nx * m; + *brows = ny * m; + *r0 = ((double) (*bcols - 1)) / 2.0; + *s0 = ((double) (*brows - 1)) / 2.0; + } else { + for (m=1; isc<0; isc++) m *= 2; /* compute power-law scale factor */ + *map_scale = base * (double) m; + *bcols = dceil( nx / (double) m); + *brows = dceil( ny / (double) m); + /* note: the following computation ensures that the lower-left corner + remains at the same location even if ny/m is a non-integer. */ + *r0 = (nx / (double) m - 1.0) / 2.0; + *s0 = (ny / (double) m - 1.0) / 2.0; + } + +} + +static int dceil(double r) +{ + int ret_val = (int)r; + if (ret_val - r > 0.0) ret_val++; + return(ret_val); +} diff --git a/src/prod/utils/src/utils.h b/src/prod/utils/src/utils.h index 900a0094..295a879b 100644 --- a/src/prod/utils/src/utils.h +++ b/src/prod/utils/src/utils.h @@ -7,6 +7,19 @@ #ifndef utils_H #define utils_H +#include + +#define UTILS_DTR ((2.0*(M_PI))/360.0) /* degrees to radians */ +#define UTILS_RTD (360.0/(2.0*(M_PI))) /* radians to degrees */ + int utils_allocate_clean_aligned_memory( void **this, size_t size ); +void utils_ease2_map_info(int iopt, int isc, int ind, + double *map_equatorial_radius_m, double *map_eccentricity, + double *e2, double *map_reference_latitude, + double *map_reference_longitude, + double *map_second_reference_latitude,double * sin_phi1, + double *cos_phi1, double *kz, + double *map_scale, int *bcols, int *brows, + double *r0, double *s0, double *epsilon); #endif // utils_H diff --git a/src/sirlib/c/Makefile b/src/sirlib/c/Makefile deleted file mode 100644 index 21eec4f4..00000000 --- a/src/sirlib/c/Makefile +++ /dev/null @@ -1,62 +0,0 @@ -# -# generic makefile for SIR format C programs/library -# written by R. Chen Feb. 2001 JPL PO.DAAC -# modified by D. Long Mar. 2014 -# -# be sure to edit include/sir3.h to reflect machine/compiler type - -TOPDIR = $(PMESDR_TOP_DIR) -LIBDIR = $(TOPDIR)/lib -INCDIR = $(TOPDIR)/include - -LIB = libcsir.a -HDRS = include/sir_geom.h - -AR = ar -CC = gcc -INSTALL = cp -f -p -MKDIR = mkdir -p -RM = rm -fr - -CONFIG_CFLAGS = -ARFLAGS = ruv -CFLAGS = -I./include -I$(INCDIR) -m64 -O3 $(CONFIG_CFLAGS) -LDFLAGS = -L. -lcsir -lm - -ifeq ($(LOCALE),ALPINEicc) - $(info Build location: $(LOCALE)) - CC = icc - AR = xiar - CFLAGS = -I./include -I$(INCDIR) $(CONFIG_CFLAGS) -O3 -ipo -g -xHost -no-prec-div # -Zp8 - ARFLAGS = crs -endif -ifeq ($(LOCALE),CUMULUSicc) - $(info Build location: $(LOCALE)) - CC = icc - AR = xiar - CFLAGS = -I./include -I$(INCDIR) $(CONFIG_CFLAGS) -O3 -ipo -g -xHost -no-prec-div # -Zp8 - ARFLAGS = crs -endif - -all : $(LIB) - -install : - $(MKDIR) $(LIBDIR) - $(INSTALL) $(LIB) $(LIBDIR) - $(INSTALL) $(HDRS) $(INCDIR) - -# create static c routines library -libcsir.a : lib/sir_geom.c - rm -f *.o - rm -f libcsir.a - $(CC) $(CFLAGS) -c lib/sir_geom.c - $(AR) $(ARFLAGS) libcsir.a *.o - $(RM) *.o - $(MKDIR) $(LIBDIR) - $(INSTALL) $(LIB) $(LIBDIR) - -clean : - $(RM) $(LIB) *.o $(LIBDIR)/$(LIB) - - - diff --git a/src/sirlib/c/include/sir_geom.h b/src/sirlib/c/include/sir_geom.h deleted file mode 100644 index c1dae62d..00000000 --- a/src/sirlib/c/include/sir_geom.h +++ /dev/null @@ -1,33 +0,0 @@ -/* - * sir_geom.h - Utilities for Calibrated Enhanced-Resolution coordinate conversions - * - * 09-Jun-2016 M. A. Hardman mhardman@nsidc.org 303-492-2969 - * Copyright (C) 2016 Regents of the University of Colorado and Brigham Young University - */ -#ifndef sir_geom_H -#define sir_geom_H - -/* Include standard headers for things like printf and __FUNCTION__ */ -#include -#include -#include - -void pixtolatlon(float x, float y, float *alon, float *alat, - int iopt, float xdeg, float ydeg, - float ascale, float bscale, float a0, float b0); -void latlon2pix(float alon, float alat, float *x, float *y, - int iopt, float xdeg, float ydeg, - float ascale, float bscale, float a0, float b0); -void ease2_map_info(int iopt, int isc, int ind, - double *map_equatorial_radius_m, double *map_eccentricity, - double *e2, double *map_reference_latitude, - double *map_reference_longitude, - double *map_second_reference_latitude,double * sin_phi1, - double *cos_phi1, double *kz, - double *map_scale, int *bcols, int *brows, - double *r0, double *s0, double *epsilon); -void f2ipix(float x, float y, int *ix, int *iy, int nsx, int nsy); - -#endif // sir_geom_H - - diff --git a/src/sirlib/c/lib/sir_geom.c b/src/sirlib/c/lib/sir_geom.c deleted file mode 100644 index 96234fd2..00000000 --- a/src/sirlib/c/lib/sir_geom.c +++ /dev/null @@ -1,952 +0,0 @@ -/* (c) 1998, 1999, 2005, 2014 BYU MERS Laboratory - - standard C BYU SIR projection transformation routines - - DGL Dec. 19, 1998 - DGL Sept. 24, 1999 + fixed polster, latlon2pix, changed args to pixtolatlon - DGL July 29, 2005 + modified and corrected EASE projection code - DGL Feb 4, 2014 + added EASE2 grid projection - DGL Feb 4, 2014 + added ease2sf - DGL Aug 2, 2014 + modified EASE2T vertical dimension from 538 to 540 - DGL Aug 12, 2014 + modified EASE2 global M25km base scale to 25025.26000 - -*/ - -#include -#include -#include - -#include "sir_geom.h" -#include "cetb.h" - -int intfix(float r) -{ - int ret_val = r; - if (ret_val - r > 0.5) ret_val--; - if (r - ret_val > 0.5) ret_val++; - return(ret_val); -} - -int dceil(double r) -{ - int ret_val = r; - if (ret_val - r > 0.0) ret_val++; - return(ret_val); -} - - -/******* standard sir format geometric transformation routines *****/ - - -void ilambert1(float*, float*, float, float, float, float, int); -void ipolster(float*, float*, float, float, float, float); -void ieasegrid(int, float*, float*, float, float, float); -void iease2grid(int, float*, float*, float, float, float, float); - -void pixtolatlon(float x, float y, float *alon, float *alat, - int iopt, float xdeg, float ydeg, - float ascale, float bscale, float a0, float b0) -{ - /* computes the lat, long position of the lower-left corner of the - x,y-th pixel. pixels indexed [1..nsx,1..nsy] */ - - float thelon, thelat; - - switch(iopt) { - case -1: - case 0: - *alon = (x-1.0)/ascale+a0; - *alat = (y-1.0)/bscale+b0; - break; - case 1: - case 2: - thelon = (x-1.0)/ascale+a0; - thelat = (y-1.0)/bscale+b0; - ilambert1(alat,alon,thelon,thelat,ydeg,xdeg,iopt); - break; - case 5: - thelon = (x-1.0)*ascale+a0; - thelat = (y-1.0)*bscale+b0; - ipolster(alon,alat,thelon,thelat,xdeg,ydeg); - break; - case 8: - case 9: - case 10: - thelon = x - 1.0 + a0; - thelat = y - 1.0 + b0; - iease2grid(iopt, alon, alat, thelon, thelat, ascale, bscale); - case 11: - case 12: - case 13: - thelon = x - 1.0 + a0; - thelat = y - 1.0 + b0; - ieasegrid(iopt, alon, alat, thelon, thelat, ascale); - break; - default: - *alon=0.0; - *alat=0.0; - } - return; -} - -void lambert1(float, float, float*, float*, float, float, int); -void polster(float, float, float*, float*, float, float); -void easegrid(int, float, float, float*, float*, float); -void ease2grid(int, float, float, float*, float*, float, float); - - -void latlon2pix(float alon, float alat, float *x, float *y, - int iopt, float xdeg, float ydeg, - float ascale, float bscale, float a0, float b0) -{ - /* computes the x,y pixel coordinates given the lat, lon position - x,y are fractional values not limited to within image */ - - static float thelon, thelat; - - switch(iopt) { - case -1: - case 0: - *x = ascale * (alon - a0)+1.0; - *y = bscale * (alat - b0)+1.0; - break; - case 1: - case 2: - lambert1(alat, alon, &thelon, &thelat, ydeg, xdeg, iopt); - *x = ascale * (thelon - a0)+1.0; - *y = bscale * (thelat - b0)+1.0; - break; - case 5: - polster(alon, alat, &thelon, &thelat, xdeg, ydeg); - *x = (thelon - a0)/ascale+1.0; - *y = (thelat - b0)/bscale+1.0; - break; - case 8: - case 9: - case 10: - ease2grid(iopt, alon, alat, &thelon, &thelat, ascale, bscale); - *x = thelon + 1.0 - a0; - *y = thelat + 1.0 - b0; - break; - case 11: - case 12: - case 13: - easegrid(iopt, alon, alat, &thelon, &thelat, ascale); - *x = thelon + 1.0 - a0; - *y = thelat + 1.0 - b0; - break; - default: - *x=0.0; - *y=0.0; - } - return; -} - -/* floating point to integer pixel location quantization routine */ - -void f2ipix(float x, float y, int *ix, int *iy, int nsx, int nsy) -{ - /* quantizes the floating point pixel location to the actual pixel value - returns a zero if location is outside of image limits - a small amount (0.002 pixels) of rounding is permitted*/ - - if (x+0.0002 >= 1.0 && x+0.0002 <= (float) (nsx+1)) - *ix = floor(x+0.0002); - else - *ix = 0; - - if (y+0.0002 >= 1.0 && y+0.0002 <= (float) (nsy+1)) - *iy = floor(y+0.0002); - else - *iy = 0; - - return; -} - - -#define DEG2RAD (double) 0.017453292 -#define RAD2DEG (double) 57.295779 - -void lambert1(float lat, float lon, float *x, float *y, - float orglat, float orglon, int iopt) -{ -/* - COMPUTES THE FORWARD TRANSFORMATION (FROM LAT/LON TO X/Y) FOR THE - LAMBERT AZIMUTHAL EQUAL-AREA PROJECTION - - SEE "MAP PROJECTIONS USED BY THE U.S. GEOLOGICAL SURVEY" - GEOLOGICAL SURVEY BULLETIN 1532, PGS 157-173 - - FOR THIS ROUTINE, A SPHERICAL EARTH IS ASSUMED FOR THE PROJECTION. THE - ERROR WILL BE SMALL FOR SMALL-SCALE MAPS. - FOR IOPT=1 A FIXED, NOMINAL EARTH RADIUS IS USED. - FOR IOPT=2 THE LOCAL RADIUS OF THE EARTH IS USED BASED ON - THE 1972 WGS ELLIPSOID MODEL (BULLETIN PG 15). - - WRITTEN BY: DGL 20 MAY 1996 (based on fortran version) - - INPUTS: - LAT (R): LATITUDE +90 TO -90 DEG WITH NORTH POSITIVE - LON (R): LONGITUDE 0 TO +360 DEG WITH EAST POSITIVE - ORGLAT (R): ORIGIN PARALLEL +90 TO -90 DEG WITH NORTH POSITIVE - ORGLON (R): CENTRAL MERIDIAN (LONGITUDE) 0 TO +360 DEG - OR -180 TO +180 WITH EAST MORE POSITIVE - IOPT (I): EARTH RADIUS OPTION - - OUTPUTS: - X,Y (R): RECTANGULAR COORDINATES IN KM -*/ - - double orglon1, lon1, x1, y1, x2, y2, c, ak, denom, eradearth; - double radearth=6378.135; /* equitorial radius of the earth */ - double f=298.26; /* 1/f based on WGS 72 model */ - double era; - - orglon1 = fmod(orglon+720.0,360.0); - lon1 = fmod(lon+720.0, 360.0); - x1 = cos(orglat*DEG2RAD); - y1 = sin(orglat*DEG2RAD); - x2 = cos(lat*DEG2RAD); - y2 = sin(lat*DEG2RAD); - c = cos((lon1-orglon1)*DEG2RAD); - -/* COMPUTE LOCAL RADIUS OF THE EARTH AT CENTER OF IMAGE */ - eradearth=6378.0; - if (iopt == 2) { /* use the local radius */ - era = (1.-1./f); - eradearth=radearth*era/sqrt(era*era*x1*x1+y1*y1); - } - - denom = 1.0 + y1*y2 + x1*x2*c; - if (denom > 0.0) - ak = sqrt(2.0/denom); - else { - fprintf(stderr,"*** Divsion error in lambert1 ***\n"); - ak = 1.0; - } - *x = (ak * x2 * sin((lon1-orglon1)*DEG2RAD) ) * eradearth; - *y = (ak * (x1*y2-y1*x2*c) ) * eradearth; - return; -} - - - -void ilambert1(float *lat, float *lon, float x, float y, - float orglat, float orglon, int iopt) -{ -/* - COMPUTES THE INVERSE TRANSFORMATION (FROM X/Y TO LAT/LON) FOR THE - LAMBERT AZIMUTHAL EQUAL-AREA PROJECTION - - SEE "MAP PROJECTIONS USED BY THE U.S. GEOLOGICAL SURVEY" - GEOLOGICAL SURVEY BULLETIN 1532, PGS 157-173 - - FOR THIS ROUTINE, A SPHERICAL EARTH IS ASSUMED FOR THE PROJECTION. THE - ERROR WILL BE SMALL FOR SMALL-SCALE MAPS. - FOR IOPT=1 A FIXED, NOMINAL EARTH RADIUS IS USED. - FOR IOPT=2 THE LOCAL RADIUS OF THE EARTH IS USED BASED ON - THE 1972 WGS ELLIPSOID MODEL (BULLETIN PG 15). - - WRITTEN BY: DGL 20 MAY 1996 (based on fortran version) - - INPUTS: - X,Y (R): RECTANGULAR COORDINATES IN KM - ORGLAT (R): ORIGIN PARALLEL +90 TO -90 DEG WITH NORTH POSITIVE - ORGLON (R): CENTRAL MERIDIAN (LONGITUDE) 0 TO +360 DEG - OR -180 TO +180 WITH EAST MORE POSITIVE - IOPT (I): EARTH RADIUS OPTION - - OUTPUTS: - LAT (R): LATITUDE +90 TO -90 DEG WITH NORTH POSITIVE - LON (R): LONGITUDE 0 TO +360 DEG WITH EAST POSITIVE -*/ - - double orglon1, rho, c, t1, t2, x1, y1, eradearth; - double radearth=6378.135; /* equitorial radius of the earth */ - double f=298.26; /* 1/f based on WGS 72 model */ - double era; - - orglon1 = fmod(orglon+720.0,360.0); - -/* COMPUTE LOCAL RADIUS OF THE EARTH AT CENTER OF IMAGE */ - eradearth=6378.0; - if (iopt == 2) { /* use the local radius */ - era = (1.-1./f); - x1=cos(orglat*DEG2RAD); - y1=sin(orglat*DEG2RAD); - eradearth=radearth*era/sqrt(era*era*x1*x1+y1*y1); - } - - x1=x/eradearth; - y1=y/eradearth; - rho=x1*x1+y1*y1; - if (rho > 0.) { - rho = sqrt(rho); - c = 2.0 * asin(rho*0.5); - *lat = asin(cos(c)*sin(orglat*DEG2RAD)+ - y1*sin(c)*cos(orglat*DEG2RAD)/rho)*RAD2DEG; - } else { - c = 0.0; - *lat = orglat; - } - - if ( !((orglat < 0.0 ? - orglat : orglat) == 90.0)) { - if (rho == 0.0) - *lon = orglon1; - else { - t1=x1*sin(c); - t2=rho*cos(orglat*DEG2RAD)*cos(c)-y1*sin(orglat*DEG2RAD)*sin(c); - *lon = orglon1 + atan2(t1,t2)*RAD2DEG; - } - } else { - if (orglat == 90.0) - *lon = orglon1 + atan2(x1,-y1)*RAD2DEG; - else - *lon = orglon1 + atan2(x1,y1)*RAD2DEG; - } - *lon = fmod(*lon+720.0, 360.0); - if (*lon > 180.0) *lon = *lon -360.0; - return; -} - - -#define dtr (double) 0.017453292 -#define abs(a) ((a) > 0.0 ? (a) : - (a)) - -double arctand(double y, double x); - -void ipolster(float *alon, float *alat, float x, float y, - float xlam, float slat) -{ -/* - COMPUTES THE INVERSE POLAR STEROGRAPHIC TRANSFORMATION FOR (X,Y) - GIVEN IN KM WITH REFERENCES LON,LAT=(XLAM,SLAT). - OUTPUT LON,LAT=ALON,ALAT - - ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES - AS RECEIVED FROM M. DRINKWATER (1994). UPDATED BY D. LONG TO - IMPROVE ACCURACY USING ITERATION WITH FORWARD TRANSFORM. -*/ - int icnt; - float xx,yy; - double r,rr,rerr,a,aa,aerr,absaerr,sn1; - double rho,e,e22,e23,cm,chi,t,x1,y1,sn,slat1; - double e2=0.006693883; - double re=6378.273; - double pi2=1.570796327; - -/* first use approximate inverse */ - - e=sqrt(e2); - e22=e2*e2; - e23=e2*e2*e2; - x1=x; - y1=y; - rho=x1*x1+y1*y1; - if (rho > 0.0) rho = sqrt(rho); - if (rho < 0.05) { - *alon = xlam; - *alat = (slat > 0.0 ? 90.0 : -90.0); - } else { - sn=1.0; - slat1=slat; - if (slat < 0.0) { - sn = -1.0; - slat1 = - slat; - } - cm=cos(slat1 * dtr)/sqrt(1.0-e2*sin(slat1 * dtr)*sin(slat1 * dtr)); - t=tan(dtr*(45.0-0.5*slat1))/ - pow((1.-e*sin(slat1*dtr))/(1.+e*sin(slat1*dtr)),e*0.5); - t=rho*t/(re*cm); - chi=pi2-2.*atan(t); - t=chi+(0.5*e2+5.*e22/24.+e23/12.0)*sin(2.*chi)+ - (7.0*e22/48.0+28.0*e23/240.0)*sin(4.*chi)+ - (7.0*e23/120.0)*sin(6.*chi); - *alat=sn*(t*90.0/pi2); - *alon=sn*atan2(sn*x1,-sn*y1)*RAD2DEG+xlam; - if (*alon < -180.0) *alon = *alon+360.0; - if (*alon > 180.0) *alon = *alon-360.0; - } - -/* using the approximate result as a starting point, iterate to improve - the accuracy of the inverse solution */ - - sn1 = 1.0; - if (slat < 0.0) sn1=-1.0; - a=arctand( (double) y, (double) x); - r=sqrt((double) (x * x+ y * y)); - icnt=0; - - label10: - polster(*alon,*alat,&xx,&yy,xlam,slat); - rr=sqrt((double) (xx*xx+yy*yy)); - rerr=sn1*(rr-r)/180.0; - aa=arctand((double) yy,(double) xx); - aerr=aa-a; - absaerr=abs(aerr); - if (absaerr > 180.0) aerr=360.0-aerr; - absaerr=abs(aerr); - if ((abs(rerr) < 0.001 && absaerr < 0.001) || icnt > 9) goto label40; - *alon = *alon +aerr; - if (abs(*alon) > 360.0) *alon=fmod(*alon,360.0); - if (*alat * slat < 0) { - rerr = rerr * (1. - sin(dtr*(*alat>0.0?*alat:- *alat))); - if ((rerr>0.0?rerr:-rerr)>2.0) - rerr=(rerr>0.0?2.0:-2.0)/(double) icnt; - } - *alat = *alat + rerr; - if (abs(*alat) > 90.0) - *alat = (*alat > 0.0 ? 90.0 : -90.0); - icnt++; - goto label10; - label40: - if (abs(*alon) > 360.0) *alon = fmod(*alon + 360.0, 360.); - return; -} - -double arctand(double y, double x) -{ - if (y == 0.0 && x == 0.0) - return(0.0); - return(atan2(y,x)*RAD2DEG); -} - -void polster( float alon, float alat, float *x1, float *y1, float xlam, float slat) -{ -/* - COMPUTES THE POLAR STEROGRAPHIC TRANSFORMATION FOR A LON,LAT - INPUT OF (ALON,ALAT) WITH REFERENCE ORIGIN LON,LAT=(XLAM,SLAT). - OUTPUT IS (X1,Y1) IN KM - - ALGORITHM IS THE SAME AS USED FOR PROCESSING ERS-1 SAR IMAGES - AS RECEIVED FROM M. DRINKWATER (1994) -*/ - double e, t, tx, ty, cm, rho, rlat, sn; - double e2=0.006693883; - double re=6378.273; - - e=sqrt(e2); - if (slat < 0.0) { - sn = -1.0; - rlat = -alat; - } else { - sn = 1.0; - rlat = alat; - } - t=pow((1.0-e*sin(dtr*rlat))/(1.0+e*sin(dtr*rlat)),e*0.5); - ty=tan(dtr*(45.0-0.5*rlat))/t; - rlat = abs(slat); - t=pow((1.0-e*sin(dtr*rlat))/(1.0+e*sin(dtr*rlat)),e*0.5); - tx=tan(dtr*(45.0-0.5*rlat))/t; - cm=cos(dtr*rlat)/sqrt(1.0-e2*sin(dtr*rlat)*sin(dtr*rlat)); - rho=re*cm*ty/tx; - *x1 = (sn*sin(dtr*(sn*alon-xlam)))*rho; - *y1 =-(sn*cos(dtr*(sn*alon-xlam)))*rho; - return; -} - - -void easegrid(int iopt, float alon, float alat, - float *thelon, float *thelat, float ascale) -{ -/* - COMPUTE THE FORWARD "EASE" GRID TRANSFORM - - GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND - THE SCALE (ASCALE) THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED - USING THE "EASE GRID" (VERSION 1.0) TRANSFORMATION GIVEN IN FORTRAN - SOURCE CODE SUPPLIED BY NSIDC - - IOPT IS EASE TYPE: IOPT=11=NORTH, IOPT=12=SOUTH, IOPT=13=CYLINDRICAL - - THE RADIUS OF THE EARTH USED IN THIS PROJECTION IS EMBEDDED INTO - ASCALE WHILE THE PIXEL DIMENSION IN KM IS EMBEDDED IN BSCALE - THE BASE VALUES ARE: RADIUS EARTH= 6371.228 KM - PIXEL DIMEN =25.067525 KM - THEN, BSCALE = BASE_PIXEL_DIMEN - ASCALE = RADIUS_EARTH/BASE_PIXEL_DIMEN - -*/ - double pi2=1.57079633; - - switch (iopt) { - case 11: /* EASE grid north */ - *thelon = ascale*sin(alon*dtr)*sin((45.0-0.5*alat)*dtr); - *thelat =-ascale*cos(alon*dtr)*sin((45.0-0.5*alat)*dtr); - break; - case 12: /* EASE grid south */ - *thelon = ascale*sin(alon*dtr)*cos((45.0-0.5*alat)*dtr); - *thelat = ascale*cos(alon*dtr)*cos((45.0-0.5*alat)*dtr); - break; - case 13: /* EASE cylindrical */ - *thelon = ascale * pi2 * alon * cos(30.0 * dtr)/90.0; - *thelat = ascale * sin(alat * dtr) / cos(30.0 * dtr); - } - return; -} - -void ieasegrid(int iopt, float *alon, float *alat, - float thelon, float thelat, float ascale) -{ -/* - COMPUTE THE INVERSE "EASE" GRID TRANSFORM - - GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND - THE SCALE (ASCALE) THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED - USING THE "EASE GRID" (VERSION 1.0) TRANSFORMATION GIVEN IN FORTRAN - SOURCE CODE SUPPLIED BY NSIDC - - IOPT IS EASE TYPE: IOPT=11=NORTH, IOPT=12=SOUTH, IOPT=13=CYLINDRICAL - - THE RADIUS OF THE EARTH USED IN THIS PROJECTION IS EMBEDDED INTO - ASCALE WHILE THE PIXEL DIMENSION IN KM IS EMBEDDED IN BSCALE - THE BASE VALUES ARE: RADIUS EARTH= 6371.228 KM - PIXEL DIMEN =25.067525 KM - THEN, BSCALE = BASE_PIXEL_DIMEN - ASCALE = RADIUS_EARTH/BASE_PIXEL_DIMEN - -*/ - double pi2=1.57079633; - double x1, y1, temp; - - x1 = thelon; - y1 = thelat; - switch (iopt) { - case 11: /* EASE grid north */ - *alon = arctand( x1, -y1); - if (abs(sin(*alon*dtr)) > abs(cos(*alon*dtr))) - temp = (x1/sin(dtr* *alon))/ ascale; - else - temp = (-y1/cos(dtr* *alon))/ ascale; - if (abs(temp) <= 1.0) - *alat = 90.0 - 2.0 * asin(temp)*RAD2DEG; - else { - *alat = -90.0 + 2.0 * asin(1./temp)*RAD2DEG; - /* - fprintf(stderr,"*** ERROR in EASE inverse sine ***\n"); - fflush(stderr); - */ - *alat = (temp > 0.0 ? 90.0 : -90.0); - } - break; - case 12: /* EASE grid south */ - *alon = arctand( x1, y1); - if (abs(cos(*alon*dtr)) > abs(sin(*alon*dtr))) - temp = (y1/cos(dtr* *alon))/ ascale; - else - temp = (x1/sin(dtr* *alon))/ ascale; - if (abs(temp) <= 1.0) - *alat = 90.0 - 2.0 * acos(temp)*RAD2DEG; - else { - /* - fprintf(stderr,"*** ERROR in EASE inverse cosine ***\n"); - fflush(stderr); - */ - *alat = 0.0; - } - break; - case 13: /* EASE cylindrical */ - *alon = ((x1/ ascale)/cos(dtr*30.0))*90.0/pi2; - temp = (y1 * cos(dtr*30.0))/ ascale; - if (abs(temp) <= 1.0) - *alat = asin(temp)*RAD2DEG; - else { - /* - fprintf(stderr,"*** ERROR in EASE inverse sine ***\n"); - fflush(stderr); - */ - *alat = (temp > 0.0 ? 90.0 : -90.0); - } - } - return; -} - -#define DTR 0.01745329241994 -#define RTD 57.29577951308232 -/* set base EASE2 map projection information */ -void ease2_map_info(int iopt, int isc, int ind, - double *map_equatorial_radius_m, double *map_eccentricity, - double *e2, double *map_reference_latitude, - double *map_reference_longitude, - double *map_second_reference_latitude,double * sin_phi1, - double *cos_phi1, double *kz, - double *map_scale, int *bcols, int *brows, - double *r0, double *s0, double *epsilon) -{ /* define key EASE2 grid information - - inputs - iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M - isc: scale factor 0..5 grid size is (basesize(ind))/2^isc - ind: base grid size index (map units per cell in m - - NSIDC .grd file for isc=0 - project type ind=0 ind=1 ind=2 - N EASE2_N25km EASE2_N36km EASE2_N24km - S EASE2_S25km EASE2_S36km EASE2_S24km - T/M EASE2_T25km EASE2_M36km EASE2_M24km - - cell size (m) for isc=0 (scale is reduced by 2^isc) - project type ind=0 ind=1 ind=2 - N 25000.0 36000.0 24000.0 - S 25000.0 36000.0 24000.0 - T/M T25025.26 M36032.220840584 M24021.4805603 - - for a given base cell size isc is related to NSIDC .grd file names - isc N .grd name S .grd name T .grd name - 0 EASE2_N25km EASE2_S25km EASE2_T25km - 1 EASE2_N12.5km EASE2_S12.5km EASE2_T12.5km - 2 EASE2_N6.25km EASE2_S6.25km EASE2_T6.25km - 3 EASE2_N3.125km EASE2_S3.125km EASE2_T3.125km - 4 EASE2_N1.5625km EASE2_S1.5625km EASE2_T1.5625km - - outputs - map_equatorial_radius_m EASE2 Earth equitorial radius (km) [WGS84] - map_eccentricity EASE2 Earth eccentricity [WGS84] - map_reference_latitude Reference latitude (deg) - map_reference_longitude Reference longitude (deg) - map_second_reference_latitude Secondary reference longitude* (deg) - sin_phi1, cos_phi1 kz EASE2 Cylin parameters* - map_scale EASE2 map projection pixel size (km) - bcols, brows, EASE2 grid size in pixels - r0, s0 EASE2 base projection size in pixels - epsilon EASE2 near-polar test factor - - *these parameters only assigned values if projection is T - - */ - - double base; - int m, nx, ny; - int region_index; - - *map_equatorial_radius_m = 6378137.0 ; /* WGS84 */ - *map_eccentricity = 0.081819190843 ; /* WGS84 */ - *e2 = *map_eccentricity * *map_eccentricity; - *map_reference_longitude = 0.0; - *epsilon = 1.e-6; - - /* map-specific parameters - these all come from cetb.h */ - switch (iopt) { - case 8: /* EASE2 grid north */ - *map_reference_latitude = 90.0; - switch(ind) { - case CETB_36KM: /* EASE2_N36km.gpd */ - region_index = CETB_36KM*CETB_NUMBER_PROJECTIONS; - base=cetb_exact_scale_m[region_index][0]; //36000.0; - nx=cetb_grid_cols[region_index][0]; //500; - ny=cetb_grid_rows[region_index][0]; //500; - break; - case CETB_24KM: /* EASE2_N24km.gpd */ - region_index = CETB_24KM*CETB_NUMBER_PROJECTIONS; - base=cetb_exact_scale_m[region_index][0]; //24000.0; - nx=cetb_grid_cols[region_index][0]; //750; - ny=cetb_grid_rows[region_index][0]; //750; - break; - default: /* EASE2_N25km.gpd */ - region_index = CETB_25KM*CETB_NUMBER_PROJECTIONS; - base=cetb_exact_scale_m[region_index][0]; //25000.0; - nx=cetb_grid_cols[region_index][0]; //720; - ny=cetb_grid_rows[region_index][0]; //720; - } - break; - case 9: /* EASE2 grid south */ - *map_reference_latitude = -90.0; - switch(ind) { - case CETB_36KM: /* EASE2_S36km.gpd */ - region_index = (CETB_36KM*CETB_NUMBER_PROJECTIONS)+1; - base=cetb_exact_scale_m[region_index][0]; //36000.0; - nx=cetb_grid_cols[region_index][0]; //500; - ny=cetb_grid_rows[region_index][0]; //500; - break; - case CETB_24KM: /* EASE2_S24km.gpd */ - region_index = (CETB_24KM*CETB_NUMBER_PROJECTIONS)+1; - base=cetb_exact_scale_m[region_index][0]; //24000.0; - nx=cetb_grid_cols[region_index][0]; //750; - ny=cetb_grid_rows[region_index][0]; //750; - break; - default: /* EASE2_S25km.gpd */ - region_index = (CETB_25KM*CETB_NUMBER_PROJECTIONS)+1; - base=cetb_exact_scale_m[region_index][0]; //25000.0; - nx=cetb_grid_cols[region_index][0]; //720; - ny=cetb_grid_rows[region_index][0]; //720; - } - break; - case 10: /* EASE2 cylindrical */ - *map_reference_latitude = 0.0; - *map_second_reference_latitude = 30.0; - *sin_phi1 = sin( DTR * *map_second_reference_latitude ); - *cos_phi1 = cos( DTR * *map_second_reference_latitude ); - *kz = *cos_phi1 / sqrt( 1.0 - *e2 * *sin_phi1 * *sin_phi1 ); - switch(ind) { - case CETB_36KM: /* EASE2_M36km.gpd */ - region_index = (CETB_36KM*CETB_NUMBER_PROJECTIONS)+2; - base=cetb_exact_scale_m[region_index][0]; //36032.220840584; - nx=cetb_grid_cols[region_index][0]; //964; - ny=cetb_grid_rows[region_index][0]; //406; - break; - case CETB_24KM: /* EASE2_M24km.gpd */ - region_index = (CETB_24KM*CETB_NUMBER_PROJECTIONS)+2; - base=cetb_exact_scale_m[region_index][0]; //24021.480560389347; - nx=cetb_grid_cols[region_index][0]; //1446; - ny=cetb_grid_rows[region_index][0]; //609; - break; - default: /* EASE2_T25km.gpd */ - region_index = (CETB_25KM*CETB_NUMBER_PROJECTIONS)+2; - base=cetb_exact_scale_m[region_index][0]; //25025.26000; - nx=cetb_grid_cols[region_index][0]; //1388; - ny=cetb_grid_rows[region_index][0]; //540; /* originally was 538 */ - } - } - - - /* grid info */ - if (isc>=0) { - for (m=1; isc>0; isc--) m *= 2; /* compute power-law scale factor */ - *map_scale = base / (double) m; - *bcols = nx * m; - *brows = ny * m; - *r0 = ((double) (*bcols - 1)) / 2.0; - *s0 = ((double) (*brows - 1)) / 2.0; - } else { - for (m=1; isc<0; isc++) m *= 2; /* compute power-law scale factor */ - *map_scale = base * (double) m; - *bcols = dceil( nx / (double) m); - *brows = dceil( ny / (double) m); - /* note: the following computation ensures that the lower-left corner - remains at the same location even if ny/m is a non-integer. */ - *r0 = (nx / (double) m - 1.0) / 2.0; - *s0 = (ny / (double) m - 1.0) / 2.0; - } - -} - -double easeconv_normalize_degrees(double b) -{ /* return -180 <= b <= 180.0 */ - while (b < -180.0) - b = b + 360.0; - while (b > 180.0) - b = b - 360.0; - return(b); -} - -void ease2sf(int iopt, float ascale, float bscale, float *fmap_scale, - int *bcols, int *brows, float *fr0, float *fs0) -{ - /* return selected EASE2 grid information (external interface) - - inputs - iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M - (ascale and bscale should be integer valued) - ascale: grid scale factor (0..5) pixel size is (bscale/2^ascale) - bscale: base grid scale index (ind=int(bscale)) - - outputs - fmap_scale EASE2 map projection pixel size (km) - bcols, brows, EASE2 grid size in pixels - fr0, fs0 EASE2 base projection size in pixels - */ - - double map_equatorial_radius_m,map_eccentricity, e2, - map_reference_latitude, map_reference_longitude, - map_second_reference_latitude, sin_phi1, cos_phi1, kz, - map_scale, r0, s0, epsilon; - - int ind = intfix(bscale); - int isc = intfix(ascale); - - /* get base EASE2 map projection parameters */ - ease2_map_info(iopt, isc, ind, &map_equatorial_radius_m, &map_eccentricity, - &e2, &map_reference_latitude, &map_reference_longitude, - &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, - &map_scale, bcols, brows, &r0, &s0, &epsilon); - *fmap_scale=map_scale; - *fr0=r0; - *fs0=s0; -} - -void ease2grid(int iopt, float alon, float alat, - float *thelon, float *thelat, float ascale, float bscale) -{ -/* - COMPUTE THE FORWARD "EASE2" GRID TRANSFORMATION - - GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND - THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED - USING THE "EASE2 GRID" (VERSION 2.0) TRANSFORMATION GIVEN IN IDL - SOURCE CODE SUPPLIED BY MJ BRODZIK - RADIUS EARTH=6378.137 KM (WGS 84) - MAP ECCENTRICITY=0.081819190843 (WGS84) - - inputs: - iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M - alon, alat: lon, lat (deg) to convert (can be outside of image) - ascale and bscale should be integer valued) - ascale: grid scale factor (0..5) pixel size is (bscale/2^ascale) - bscale: base grid scale index (ind=int(bscale)) - - outputs: - thelon: X coordinate in pixels (can be outside of image) - thelat: Y coordinate in pixels (can be outside of image) - -*/ - double map_equatorial_radius_m,map_eccentricity, e2, - map_reference_latitude, map_reference_longitude, - map_second_reference_latitude, sin_phi1, cos_phi1, kz, - map_scale, r0, s0, epsilon; - int bcols, brows; - - int ind = intfix(bscale); - int isc = intfix(ascale); - double dlon = alon; - double phi = DTR * alat; - double lam = dlon; - - double sin_phi, q, qp, rho, x, y; - - /* get base EASE2 map projection parameters */ - ease2_map_info(iopt, isc, ind, &map_equatorial_radius_m, &map_eccentricity, - &e2, &map_reference_latitude, &map_reference_longitude, - &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, - &map_scale, &bcols, &brows, &r0, &s0, &epsilon); - - dlon = dlon - map_reference_longitude; - dlon = easeconv_normalize_degrees( dlon ); - lam = DTR * dlon; - - sin_phi=sin(phi); - q = ( 1.0 - e2 ) * ( ( sin_phi / ( 1.0 - e2 * sin_phi * sin_phi ) ) - - ( 1.0 / ( 2.0 * map_eccentricity ) ) - * log( ( 1.0 - map_eccentricity * sin_phi ) - / ( 1.0 + map_eccentricity * sin_phi ) ) ); - - switch (iopt) { - case 8: /* EASE2 grid north */ - qp = 1.0 - ( ( 1.0 - e2 ) / ( 2.0 * map_eccentricity ) - * log( ( 1.0 - map_eccentricity ) - / ( 1.0 + map_eccentricity ) ) ); - if ( abs( qp - q ) < epsilon ) - rho = 0.0; - else - rho = map_equatorial_radius_m * sqrt( qp - q ); - x = rho * sin( lam ); - y = -rho * cos( lam ); - break; - case 9: /* EASE2 grid south */ - qp = 1.0 - ( ( 1.0 - e2 ) / ( 2.0 * map_eccentricity ) - * log( ( 1.0 - map_eccentricity ) - / ( 1.0 + map_eccentricity ) ) ); - if ( abs( qp + q ) < epsilon ) - rho = 0.0; - else - rho = map_equatorial_radius_m * sqrt( qp + q ); - x = rho * sin( lam ); - y = rho * cos( lam ); - break; - case 10: /* EASE2 cylindrical */ - x = map_equatorial_radius_m * kz * lam; - y = ( map_equatorial_radius_m * q ) / ( 2.0 * kz ); - break; - default: - fprintf(stderr,"*** invalid EASE2 projection specificaion %d in ease2grid\n",iopt); - break; - } - - *thelon = (float) (r0 + ( x / map_scale ) + 0.5); - *thelat = (float) (s0 + ( y / map_scale ) + 0.5); /* 0 at bottom (BYU SIR files) */ - - return; -} - -void iease2grid(int iopt, float *alon, float *alat, - float thelon, float thelat, float ascale, float bscale) -{ -/* - COMPUTE THE INVERSE "EASE2" GRID TRANSFORM - - GIVEN THE IMAGE TRANSFORMATION COORDINATES (THELON,THELAT) AND - THE CORRESPONDING LON,LAT (ALON,ALAT) IS COMPUTED - USING THE "EASE GRID" (VERSION 2.0) TRANSFORMATION GIVEN IN IDL - SOURCE CODE SUPPLIED BY MJ BRODZIK - - inputs: - iopt: projection type 8=EASE2 N, 9-EASE2 S, 10=EASE2 T/M - thelon: X coordinate in pixels (can be outside of image) - thelat: Y coordinate in pixels (can be outside of image) - ascale and bscale should be integer valued) - ascale: grid scale factor (0..5) pixel size is (bscale/2^ascale) - bscale: base grid scale index (ind=int(bscale)) - - outputs: - alon, alat: lon, lat location in deg (can be outside of image) -*/ - double map_equatorial_radius_m,map_eccentricity, e2, - map_reference_latitude, map_reference_longitude, - map_second_reference_latitude, sin_phi1, cos_phi1, kz, - map_scale, r0, s0, epsilon; - int bcols, brows; - - int ind = intfix(bscale); - int isc = intfix(ascale); - - double lam, arg, phi, beta, qp, rho2, x, y, e4, e6; - - /* get base EASE2 map projection parameters */ - ease2_map_info(iopt, isc, ind, &map_equatorial_radius_m, &map_eccentricity, - &e2, &map_reference_latitude, &map_reference_longitude, - &map_second_reference_latitude, &sin_phi1, &cos_phi1, &kz, - &map_scale, &bcols, &brows, &r0, &s0, &epsilon); - e4 = e2 * e2; - e6 = e4 * e2; - - /* qp is the function q evaluated at phi = 90.0 deg */ - qp = ( 1.0 - e2 ) * ( ( 1.0 / ( 1.0 - e2 ) ) - - ( 1.0 / ( 2.0 * map_eccentricity ) ) - * log( ( 1.0 - map_eccentricity ) - / ( 1.0 + map_eccentricity ) ) ); - - x = ((double) thelon - r0 - 0.5) * map_scale; - //y = (s0 - (double) thelat + 0.5) * map_scale; /* 0 at top (NSIDC) */ - y = ((double) thelat - 0.5 - s0) * map_scale; /* 0 at bottom (BYU SIR files) */ - - switch (iopt) { - case 8: /* EASE2 grid north */ - rho2 = ( x * x ) + ( y * y ); - arg=1.0 - ( rho2 / ( map_equatorial_radius_m * map_equatorial_radius_m * qp ) ); - if (arg > 1.0) arg=1.0; - if (arg < -1.0) arg=-1.0; - beta = asin( arg ); - lam = atan2( x, -y ); - break; - case 9: /* EASE2 grid south */ - rho2 = ( x * x ) + ( y * y ); - arg = 1.0 - ( rho2 / ( map_equatorial_radius_m * map_equatorial_radius_m * qp ) ); - if (arg > 1.0) arg=1.0; - if (arg < -1.0) arg=-1.0; - beta = -asin( arg ); - lam = atan2( x, y ); - break; - case 10: /* EASE2 cylindrical */ - arg = 2.0 * y * kz / ( map_equatorial_radius_m * qp ); - if (arg > 1.0) arg=1.0; - if (arg < -1.0) arg=-1.0; - beta = asin( arg ); - lam = x / ( map_equatorial_radius_m * kz ); - break; - default: - fprintf(stderr,"*** invalid EASE2 projection specification %d in iease2grid\n",iopt); - break; - } - - phi = beta - + ( ( ( e2 / 3.0 ) + ( ( 31.0 / 180.0 ) * e4 ) - + ( ( 517.0 / 5040.0 ) * e6 ) ) * sin( 2.0 * beta ) ) - + ( ( ( ( 23.0 / 360.0 ) * e4) - + ( ( 251.0 / 3780.0 ) * e6 ) ) * sin( 4.0 * beta ) ) - + ( ( ( 761.0 / 45360.0 ) * e6 ) * sin( 6.0 * beta ) ); - - *alat = (float) (RTD * phi); - *alon = (float) (easeconv_normalize_degrees( map_reference_longitude + ( RTD*lam ) ) ); - - return; -} - diff --git a/src/sirlib/c/readme_c.txt b/src/sirlib/c/readme_c.txt deleted file mode 100644 index 2cb1c7a0..00000000 --- a/src/sirlib/c/readme_c.txt +++ /dev/null @@ -1,60 +0,0 @@ -C routines for the BYU-MERS "SIR" image format - - -The BYU-MERS "sir" image format was developed by the Brigham Young -University (BYU) Microwave Earth Remote Sensing (MERS) research group -to store images of the earth along with the information required to -earth-locate the image pixels. - -Note that these routines use the C-standard pixel address from 0 to N-1 -rather than the SIR-standard pixel address from 1 to N. Note, however, that -(1,1) is in the lower left of the image of a SIR File. For an input location -of pixel (1,1) the Lat,Lon values returned by pixtolatlon correspond to -the location of the lower-left corner of the (1,1)th pixel. - -This directory contains several programs which illustrate reading and -writing BYU SIR files. The programs include utilities to convert the -SIR files into other file formats. - -Library routines are in lib/ and include sir_io.c which has the read/ -write routines, sir_geom.c which contains the geometry routines, and -sir_ez.c which is an easy interface to the low-level routines in -the other two files. The library interfaces are documented in the -include files include/sir3.h (for sir_io.c and sir_geom.c) and -include/sir_ez.h for the sir_ez.c routines. - -Programs illustrating of how to call the various routines and to convert -SIR files to various other image types are in this directory. Given -an input SIR file, the program csir_locmap creates SIR images of the -latitude and longitudue values of the SIR file. sir_dump and sir_dump_small -create an ASCII file of the SIR image values and locations. - -The tools/ directory contains additional example programs and utilities to -land mask SIR files (sirmask), compute the difference of two SIR files -(sirdiff) extract subimages to another SIR file (sir_extractregion), modify -the heading of a SIR file (sir_modhead), and combine and remap SIR images into -different SIR projections (sir_remapper and sir_remapper2). Generally, it -is recommended that header info be initialized from a previously read in \ -sir file and modified as necessary. - -If you have any questions, please contact me. - - -============================================================================== - -Dr. David G. Long long@ee.byu.edu -Director, BYU Center for Remote Sensing www.cers.byu.edu -Professor, Electrical and Computer Eng. Dept. www.ee.byu.edu -Brigham Young University www.byu.edu -459 Clyde Building voice: 801-422-4383 -Provo, Utah 84602 fax: 801-422-0201 - -Scatterometer Climate Record Pathfinder: www.scp.byu.edu - -============================================================================== - -Code is SIR header version 3.0 compliant. - -Last revised: 26 Feb 2014 DGL -(c) 1999, 2000, 2002, 2014 BYU MERS -