[McIDAS #UQI-731403]: Remap2

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);
*/
}