Hi Randy and Mary Ellen, Previously I said: > I note a number of differences in the SSEC REMAP2 code and the REMAP2_TASC > code; some look like they could be important, and some not (REMAP2_TASC > is obviously based on an earlier version of REMAP2). It seems to me that > it would be a good idea to try and fit your 'cloudp' code into a modified > copy of the SSEC REMAP2 code. The reason I say this is that the SSEC code > has been modified much more recently than the code on which REMAP2_TASC > was based. I just went finished revising the remap2_tasc.pgm code to match the SSEC XRD v2012.1 code everywhere possible. The major difference in the two routines should now be: - remap2_tasc.pgm supports user-specified LAT,LON bounds and uses those values to create the RECT navigation in the output image - remap2_tasc.pgm calls the cloudp routine to modify band values I have attached the modified remap2_tasc.pgm code to this email so you can give it a test in your environment. I also modified cloudp_tasc.c to set the variable 'five' to '5' unless the satellite sensor source is >= 180. This should allow the routine to be run without modification for both GOES-11 (and previous) and GOES-12 (and later) imagery without the need to modify the cloudp_tasc.c code by hand. I have attached the modified cloudp_tasc.c code. Comment: - I can remap AREA9190 and AREA9191 with no problems, but the values in the remapped bands are all zero for bands 1, 2, 4 and 5/6 and 255 for band 3. Do you see this? Cheers, Tom -- **************************************************************************** Unidata User Support UCAR Unidata Program (303) 497-8642 P.O. Box 3000 address@hidden Boulder, CO 80307 ---------------------------------------------------------------------------- Unidata HomePage http://www.unidata.ucar.edu **************************************************************************** Ticket Details =================== Ticket ID: UQI-731403 Department: Support McIDAS Priority: Normal Status: Closed
Attachment:
remap2_tasc.pgm
Description: Portable graymap
#include <math.h>
#include <stdio.h>
/* Function to calculate Reflective component of the 3.9 mm */
/* channel during daytime correcting for Solar Zenith Angle */
/* Inputs are 3.9mm channel radiance and solar zenith angle */
/* C1 (mW/(m2.ster.cm-1), C2 (cm K)are Planck function constants */
/* R = radius of sun (m) orbit (m) */
/* r = radius of earth's obrit (m) */
/* Nsun = radiance of blackbody temperature of solar photosphere */
/* S = correction for solar zenith angle */
/* Sza = solar zenith angle */
/* PlanckRadiance = emitted component for a blackbody at temp (K) */
/* found by computing the planck function. Input is Channel 4 temp */
/* Reflec = Reflected component of 3.9 mm channel radiance */
/* N = 3.9mm radiance */
/* Lambda = 3.9 microns WN (wavenumber) = 2564.10 cm-1 */
/* NOTE User must use a strech table for each channel in order to view
imagery. Use SU INI "NAME" GVAR RAW
SU MAKE "NAME" lower value upper value 0 255 */
/* MEL 11/20/02 change sza def for ref product from <76deg to <88deg */
extern float ftime_(int*);
double MinReflec=100.0, MaxReflec=0.0; /* <<<<< UPC mod 20030610 >>>>> */
int cnt=0;
#if 0 /* <<<<< UPC mod 20030621 >>>>> */
#include "sza_tasc.c"
#endif
long julday(int y,int m,int d)
{long jul;
int ja,jy=y,jm;
if(m>2) jm=m+1;
else
{--jy;
jm=m+13;
}
jul=(long)(floor(365.25*jy)+floor(30.6001*jm)+d+1720995);
ja=(int)(0.01*jy);
jul+=2-ja+(int)(0.25*ja);
return jul;
}
double vdot(double *x, double *y)
{int i;
double sum;
sum=0.0;
for(i=0;i<3;i++)sum+=x[i]*y[i];
return sum;
}
float sun_el(int y, int m, int d, float h, float lat, float lon)
{long jd;
int i;
double rlat,rlon,coslat,costheta,ddmag,ecc,denom;
double rjd,gha,tn,tu,g,lambda,epsilon,smlon,sr,sun_t[3],surf[3],dd[3];
rlat=(double)lat*M_PI/180.0;
rlon=(double)lon*M_PI/180.0;
jd=julday(y,m,d);
rjd=(double)jd+(double)(h-12)/24.0;
tn=rjd-2451545.0;
tu=(double)(jd-2451545)/36525;
gha=24110.54841+(8640184.812866+(0.093104-6.2e-6*tu)*tu)*tu;
gha=360.0*gha/86400.0+(double)h*15.0;
while(gha<0.0)gha+=360.0;
gha*=M_PI/180.0;
smlon=280.466+0.9856474*tn;
g=357.528+0.9856003*tn;
g*=M_PI/180.0;
epsilon=23.440-0.0000004*tn;
epsilon*=M_PI/180.0;
sr=1.00014-0.01671*cos(g)-0.00014*cos(2.0*g);
lambda=smlon+1.915*sin(g)+0.020*sin(2.0*g);
lambda*=M_PI/180.0;
sr*=23454.791;
sun_t[0]=sr*cos(lambda);
sun_t[1]=sr*cos(epsilon)*sin(lambda);
sun_t[2]=sr*sin(epsilon)*sin(lambda);
ecc=0.0818173;
denom=sqrt(1.0-ecc*ecc*sin(rlat)*sin(rlat));
coslat=cos(rlat);
surf[0]=coslat*cos(rlon+gha)/denom;
surf[1]=coslat*sin(rlon+gha)/denom;
surf[2]=(1.0-ecc*ecc)*sin(rlat)/denom;
for(i=0;i<3;i++)dd[i]=sun_t[i]-surf[i];
ddmag=sqrt(vdot(dd,dd));
for(i=0;i<3;i++)dd[i]/=ddmag;
costheta=vdot(dd,surf);
return (float) acos(costheta)*180.0/M_PI;
}
void cloudp_(short *data, int *offset, int *line, int *elem, int *areadir)
{
double C1 = 1.1910439e-5, C2 = 1.438769, R = 6.9595e8, r = 1.4956e11,
lambda = 3.9;
double WN = 2564.10;
double Nsun = 2.258e5;
double Solar2Earth = 2.16533e-5;
double S, ch4_radiance, ch2_radiance,N, Reflec, Reflec1, Sza, ch4_temp,
ch3_temp;
double ch1_albedo, ch2_temp, ch5_temp,deltaT;
double dum,MaxSza, MinSza;
short ch1_buf[2], ch2_buf[2], ch4_buf[2], ch5_buf[2], temp, temp1;
int status, one, two, four, five, JDay, Year,day,mo;
int ch1=0, ch2=1, ch3=2, ch4=3, ch5=4;
float satline, satelem, dummy1, dummy2, lat, lon,SzaJerry;
char option[4];
float angles[3],Hr;
float input[4];
int dummy[40];
static int count=0;
one = 1;
two = 2;
four= 4;
/*
** <<<<< UPC mod 20120725 - for goes-12/13/14/15 change five=5 to five=6 >>>>>
** NB: - this code assumes GOES data!
*/
if ( areadir[2] >= 180 ) {
five = 6;
} else {
five= 5;
}
if ( count == 0 ) {
Mcdprintf("ss = %d, five = %d\n", areadir[2], five);
count = 1;
}
ch1_buf[0] = data[(*offset)/2 + 0];
ch2_buf[0] = data[(*offset)/2 + 1];
ch4_buf[0] = data[(*offset)/2 + 3];
ch5_buf[0] = data[(*offset)/2 + 4];
/* DUMMY array send to kb1cal and kb3cal inorder to run the new calibration in
Mcidas7.5 rja 12/8/98 */
/* <<<<< UPC mod 20120110 - pass 'dummy' as first parameter to kb2cal >>>>> */
status = kb2cal_(dummy,areadir, &one, &one, ch1_buf);
status = kb1cal_(dummy,areadir, &one, &two, ch2_buf);
status = kb1cal_(dummy,areadir, &one, &four, ch4_buf);
status = kb3cal_(dummy,areadir, &one, &five, ch5_buf);
ch1_albedo = ((double) ch1_buf[0])/10.;
ch2_temp = ((double) ch2_buf[0])/10.;
ch4_temp = ((double) ch4_buf[0])/10.;
ch5_temp = ((double) ch5_buf[0])/10.;
/*
Mcdprintf("line %d, elem %d, ch2 temp %lf, ch4 temp %lf\n", *line, *elem,
ch2_temp, ch4_temp);
*/
satline = (float) *line;
satelem = (float) *elem;
status = nv1sae_(&satline, &satelem, &dummy1, &lat, &lon, &dummy2);
strncpy(option, "ANG ", 4);
input[0] = (float) areadir[3];
input[1] = ftime_(&areadir[4]);
input[2] = lat;
input[3] = -lon;
/*
Mcprintf("%f line %d elem %d ",input[1],*line,*elem);
*/
/* calculate angles from mcidas routines */
status = nv1opt_(option, input, angles);
/*
Mcdprintf("angles %lf %lf %lf\n", angles[0], angles[1], angles[2]);
*/
Sza = angles[1];
/*
Include Jerry's routine to calculate solar zenith angle given time/day lat/lon
Year = (int) input[0]/1000 + 1900;
JDay = (int) input[0]%1000;
day=0;
for(mo=1; mo < 13; mo++)
{
day+= 30 + ((mo+(mo>>3))&1) - (1+((Year&3)!=0))*(mo==2);
if (day >= JDay) break;
}
day = (30 + ((mo+(mo>>3))&1) - (1+((Year&3)!=0))*(mo==2)) - (day-JDay);
Hr = input[1];
SzaJerry = sun_el(Year,mo,day,Hr,lat,-lon);
Mcprintf(" Sza: %lf %lf %lf \n",Sza,SzaJerry,Sza-SzaJerry);
*/
/* Calculate */
if ( (Sza >= 0.) && (Sza < 88.) )
{
/* Daytime: Calculate Reflectivity */
S = Nsun * Solar2Earth * (double)cos(M_PI/180.0*Sza);
if ( (ch4_temp != 0.) && (ch2_temp != 0.0) )
{
double num,den2,den4;
/*
Mcprintf(" C1= %lf C2=%lf WN=%lf ch2_tmp=%lf ch4_tmp=%lf
\n",C1,C2,WN,ch2_temp,ch4_temp);
*/
num=C1*pow(WN,3.0);
den4=exp(C2*WN/ch4_temp)-1.;
den2=exp(C2*WN/ch2_temp)-1.;
ch4_radiance = num/den4;
ch2_radiance = num/den2;
/*
Mcprintf(" num=%lf den2=%lf ch2_rad= %lf den4=%lf ch4_rad=%lf\n",num, den2,
ch2_radiance,den4,ch4_radiance);
*/
/* correct for solar zenith angle */
/* MEL 10.08.99 insure that Reflec calculation denominator can't be equal to 0
*/
if (S != ch4_radiance ){
Reflec = (ch2_radiance-ch4_radiance) / (S - ch4_radiance);
}else {
Reflec = MinReflec;
}
if (Reflec < MinReflec) MinReflec = Reflec;
if (Reflec > MaxReflec) MaxReflec = Reflec;
temp = 245.0*(Reflec + 0.1)/1.7;
if (temp < 10) temp=10;
if (temp > 255) temp=255;
data[(*offset)/2 + ch3] = temp; /* Set Channel 3 daytime reflec */
} else {
data[(*offset)/2 + ch3] = 255;
}
/*
Mcprintf("lat=%lf lon=%lf ch2:%5.4lf %5.2lf ch4:%5.4lf %5.2lf Sza:%5.2lf
S=%8.4lf Reflec: %lf Scaled:%d\n",lat,-lon,ch2_radiance,ch2_temp,ch4_radiance,
ch4_temp,Sza,S, Reflec, temp);
*/
}
else if (Sza >= 88.0)
{
/* Nighttime: Calculate Channel 4 Temp - Channel 2 Temp */
deltaT = ch4_temp - ch2_temp;
data[(*offset)/2 + ch3] = 128 + (short) (deltaT*10.); /* Set Ch3 to Nighttime
FOG */
temp=data[(*offset)/2 + ch3];
/*
Mcprintf("lat=%lf lon=%lf ch2: %lf ch4: %lf Sza:%5.2lf deltaT: %lf Scaled:%d
\n",lat,-lon,ch2_temp, ch4_temp,Sza, deltaT, temp);
*/
}
else
{
/* TERMINATER REGION SET CH 2 and 3 to 0.0 */
data[(*offset)/2 + ch3] = 0;
temp=0;
}
data[(*offset)/2 + ch1] = (short) (ch1_albedo*10.); /* channel 1 albedo */
data[(*offset)/2 + ch2] = (short) (ch2_temp*10.); /* Channel 2 */
data[(*offset)/2 + ch4] = (short) (ch4_temp*10.); /* Channel 4 */
data[(*offset)/2 + ch5] = (short) (ch5_temp*10.); /* channel 5 */
/*
Mcprintf("lat=%5.2lf lon=%5.2lf %5.2lf %5.2lf %d %5.2lf %5.2lf
\n",lat,-lon,ch1_albedo,ch2_temp,temp,ch4_temp,ch5_temp);
*/
}