[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
20000808: Converting AVHRR data to AREA form with conlev1b.c
- Subject: 20000808: Converting AVHRR data to AREA form with conlev1b.c
- Date: Tue, 08 Aug 2000 12:07:02 -0600
>From: "Larry D. Oolman" <address@hidden>
>Organization: University of Wyoming
>Keywords: 200008081721.e78HLUT05601 McIDAS-X conlev1b.c
Larry,
> Is the convert_lev1b program from the mcidas XRD package available?
>I have some AVHRR data I would like to look at.
I am putting this into my 7.7 distribution. To get you going quickly,
I have included the source for it below (conlev1b.c). Right off of the
top of my head, I can't remember if this needs any supporting routines,
so let me know if you have problems building it with McIDAS 7.6.
Here is some useful commentary from John Benson of SSEC to another user
regarding use of this conversion routine. Note the final comment:
From: John Benson <address@hidden>
Subject: Re: Converting POES Level 1b to McIDAS Format
Date: Mon, 10 Apr 2000 10:59:25 -0500 (CDT)
Here's a quick-check.
The program knows that the 4th byte of the header record is an ascii
blank (hex 20). It checks the fourth byte of the file. If it isn't
right, it assumes the first 512 bytes of the file are an SAA preamble.
So it tabs 512 bytes into the file and tries again. If it still can't
find a blank, it gives up.
So take a look at your raw data. The first four bytes should probably
be N S S <blank>. Either starting at the very beginning of the file,
or starting at byte 512. If they aren't, we don't know what you've
got.
The next thing to check is the record length. If it's the right
format, each record (including the header) should all be the same
length, which is 22,528 bytes.
--address@hidden
From: John Benson <address@hidden>
Subject: Re: Converting POES Level 1b to McIDAS Format
To: address@hidden
Date: Mon, 10 Apr 2000 13:30:03 -0500 (CDT)
I find a proper header record at data + 512 bytes, which indicates that
the data is GAC. That means the record length should be 5632 bytes.
But when I look at byte 6144 (5632 + 512) I do not see reasonable
data. The date is somewhere in the sixth century C.E.
So I'm not sure exactly what you've got.
--address@hidden
11-APR-2000 16:01 deew
Told Greg that there is a command in XRD, CONVERT_LEV1B, that
works on NOAA-15 level 1B data.
Greg was unable to get it to work.
Contacted Greg to get the specifics of his request. We did some testing
here and found that you MUST REQUEST 16 bit data and all bands. If you
request 8 or 10 bit and fewer bands the program in XRD will not make an
usable area.
------------------------------ conlev1b.c ------------------------------------
/*
* Copyright(c) 1997, Space Science and Engineering Center, UW-Madison
* Refer to "McIDAS Software Acquisition and Distribution Policies"
* in the file mcidas/data/license.txt
*/
/**** $Id: convert_lev1b.c,v 1.2 1999/05/14 16:07:39 rickk Tst $ ****/
/*
*? CONLEV1B - Converts SAA Level 1B data to McIDAS AREA format
*? CONLEV1B source_file area_number
*? Parameters:
*? source_file | file name containing level 1B AVHRR data
*? area_number | McIDAS area number to store AVHRR data
*? Remarks:
*? Note that Level 1B data is written in different formats.
*? CONLEV1B works on data written by SAA.
*? ----------
*/
#include <sys/types.h>
#include <string.h>
#include <stdio.h>
#include <unistd.h>
#include <fcntl.h>
#include <stdlib.h>
#include <errno.h>
int i4(unsigned char *);
/*
* AREA header offsets
*/
#define DATA (256+NAVSIZE+CALSIZE) /* area byte of first data line */
#define NAVSIZE 512
#define CALSIZE 512 /* */
/*
* area line offsets
*/
#define DOC 0 /* number of bytes per line for doc */
#define CAL 212 /* number of bytes per line for calibration */
/*
* the following is an extraction from the NOAA/NESDIS Level 1b
* AVHRR/3 Format document
*
*
*/
#define REC 5632 /* length of 1B record */
#define LACLEN 22528
#define RECLEN (GAC ? REC : LACLEN)
#define PIXELS (GAC ? 408 : 2048)
#define CHANLS 5
#define element(x) (1264+10*(x))
/*
* this batch of words are found in the input record
*/
#define SCAN(x) (x[1]+256*x[0])
#define SOUTH(x) (x[12]&0x80)
#define CH3(x) ((x[13]&3)==0) /* true means 3b */
#define STIM(x) ((int)(x[11]+256*(x[10]+256*(x[9]+256*x[8]))))
#define Q1(x) ((unsigned int)(x[27]+256*(x[26]+256*(x[25]+256*x[24]))))
#define Q2(x) ((unsigned int)(x[31]+256*(x[30]+256*(x[29]+256*x[28]))))
#define BAD(x) (Q1(x)&0x80000000)
#define LOSS(x) (Q1(x)&0x20000000)
#define NOCAL(x) (Q2(x)&0x0000a000)
/*
* all these words are found in the header record
* which is read into a global array when thigs get started
* so they are always available
*/
#define SID (head[73]+256*head[72])
#define TYP (head[77]+256*head[76])
#define GAC (TYP==2)
#define AVHR (TYP<4)
#define YY (head[85]+256*head[84])
#define DAY (head[87]+256*head[86])
#define MS ((int)(head[91]+256*(head[90]+256*(head[89]+256*head[88]))))
#define SS ((MS/1000)%60)
#define MM ((MS/60000)%60)
#define HH ((MS/3600000))
/* orbit epoch constants */
#define EY ((head[349]+256*head[348])%100)
#define ED (head[351]+256*head[350])
#define ET i4(&head[352])
#define ES ((ET/1000)%60)
#define EM ((ET/60000)%60)
#define EH ((ET/3600000))
/* orbital elements at epoch time */
#define SMA (i4(&head[356]))
#define ECC (i4(&head[360]))
#define INC (i4(&head[364]))
#define PER (i4(&head[368]))
#define RAS (i4(&head[372]))
#define MAN (i4(&head[376]))
#define CAL1 48
#define CAL2 108
#define CAL3a 168
#define CAL3b 228
#define CAL4 252
#define CAL5 276
#define RADCOEF 256
int area(int);
void area_lines(int,int);
void do_line(unsigned char *,unsigned short*);
void write_line(int,int,unsigned char *,int);
void load_cal(unsigned char *,char *);
void load_cal_hdr(void);
void avhr_nav(void);
extern char *Mcpathname(char *);
int noaa_number(int);
int lines=-1;
int south=0;
int ch3present=0;
int ch6present=0;
int start_time=0;
int start_scan;
unsigned char head[REC];
static
unsigned char input[LACLEN];
static
unsigned char output[CAL+DOC+2048*CHANLS*2];
int
main(int argc,char **argv)
{
int fdin,fdout;
char *arg1,*arg2;
Mcinit(argc,argv);
Mccmdstr(0,1," ",&arg1);
fdin = open(Mcpathname(arg1),O_RDONLY);
if(fdin<0)
{
perror(arg1);
return 3;
}
read(fdin,head,REC);
if(head[3] != ' ')
{
lseek(fdin,512,SEEK_SET);
read(fdin,head,REC);
if(head[3] != ' ')
{
printf("file doesn't seem to be level-1b format\n");
return 1;
}
}
if( ! AVHR)
{
printf("file is not AVHR data\n");
return 2;
}
if(!GAC)lseek(fdin,LACLEN-REC,SEEK_CUR);
if(!GAC) printf("LAC/HRPT\n");
Mccmdstr(0,2,"1",&arg2);
fdout = area(atoi(arg2));
load_cal_hdr();
while(1)
{
int threeor6;
int c = read(fdin,input,RECLEN);
if(c!=RECLEN) break;
threeor6=CH3(input);
if(lines<0)
{
/*
* overlay timestamps in header
* with those from first line
*/
memcpy(&head[84],&input[2],4);
memcpy(&head[88],&input[8],4);
start_time=STIM(input);
if(SOUTH(input)) south=1;
start_scan=SCAN(input);
}
if(BAD(input)!=0) { continue; }
if(NOCAL(input)!=0) { continue; }
if(SCAN(input)<start_scan) { continue ;}
lines = SCAN(input)-start_scan;
load_cal(input,(char *)output);
do_line(input,(unsigned short *)&output[CAL+DOC]);
write_line(fdout,lines,output,threeor6);
if(threeor6)
{
ch3present=1;
}
else
{
ch6present=1;
}
}
area_lines(fdout,lines);
close(fdout);
printf("AVHR processing done\n");
return 0;
}
static int dir[64];
void hirs_nav(void);
int nav[NAVSIZE/4];
int cal[CALSIZE/4];
/*
* modify an existing area
* by changing the number of lines
* and not incidently, adding navigation
*
* it is done last, right before the area is closed
*/
void area_lines(int fd, int lines)
{
(void)lseek(fd,(off_t)0,SEEK_SET);
(void)read(fd,dir,sizeof(dir));
if(ch6present) dir[18]|=0x20;
if(ch3present) dir[18]|=0x04;
dir[3]=DAY+1000*(YY-1900);
dir[4]=10000*HH+100*MM+SS;
dir[8]=lines;
(void)lseek(fd,(off_t)0,SEEK_SET);
(void)write(fd,dir,sizeof(dir));
avhr_nav();
write(fd,nav,NAVSIZE);
write(fd,cal,CALSIZE);
}
/*
* return fd of an area you just made, positioned ready for writing
*
*/
int area(int number)
{
char name[32];
int fd;
sprintf(name,"AREA%4.4d",number);
fd=open(Mcpathname(name),O_CREAT | O_TRUNC | O_RDWR,0644);
memset(dir,0,sizeof(dir));
memset(nav,0,sizeof(nav));
memset(cal,0,sizeof(cal));
dir[1]=4;
dir[2]=50+noaa_number(SID);
dir[3]=DAY+1000*(YY-1900);
dir[4]=10000*HH+100*MM+SS;
dir[5]=1;
dir[6]=1;
dir[7]=1;
dir[8]=0; /* lines, will change later */
dir[9]=PIXELS;
dir[10]=2;
dir[11]=1;
dir[12]=1;
if(GAC)
{
dir[11]=3;
dir[12]=5;
}
dir[13]=5; /* bands */
dir[16]=dir[3];
dir[17]=dir[4];
dir[18]=0x1b; /* 1,2,4,5 */
strcpy((char *)&dir[24],"AVHR/3 from level 1b ");
dir[33]=DATA; /* offset to data */
dir[34]=sizeof(dir); /* offset to nav */
dir[53]=2; /* sample/average = neither */
dir[62]=256+NAVSIZE; /* offset to cal */
dir[48]=DOC; /* doc */
dir[49]=CAL; /* cal */
dir[50]=8; /* lev */
dir[35]=-1; /* val */
strcpy((char *)&dir[51],"AVH3");
strcpy((char *)&dir[52],"RAW ");
dir[14]=dir[48]+dir[49]+dir[50]+4;
write(fd,dir,sizeof(dir));
write(fd,nav,NAVSIZE); /*cal & nav space */
write(fd,cal,CALSIZE);
return fd;
}
static
void re_order(unsigned char *inp, unsigned short *oup)
{
int i;
for(i=0;i<5;i++)
{
short s = inp[1]+256*inp[0];
inp+=2;
oup[i] = s;
}
}
void do_line(unsigned char *buf, unsigned short *out)
{
int i;
for(i=0;i<PIXELS;i++)
{
re_order(&buf[element(i)], out);
out += 5;
}
}
void write_line(int fd, int line, unsigned char *p,int ch3)
{
static unsigned char
order[]={1,2,3,4,5,0,0,0};
int length = (DOC+4+sizeof(order)+CAL+PIXELS*CHANLS*2);
int valcode = -1;
if(ch3)
{
order[2]=3;
}
else
{
order[2]=6;
}
lseek(fd,DATA+line*(length),SEEK_SET);
write(fd,&valcode,4);
if(DOC)write(fd,&p[0],DOC);
if(CAL)write(fd,&p[DOC],CAL);
write(fd,order,sizeof(order));
write(fd,&p[DOC+CAL],CHANLS*2*PIXELS);
}
static int days[] = { 0, 31, 28, 31, 30, 31, 30,
31, 31, 30, 31, 30, 31};
void
avhr_nav(void)
{
memset(nav,0,NAVSIZE);
if(EY%4==0) days[2]=29;
strcpy((char *)nav,"TIRO");
nav[1]=DAY + 1000*(YY%100) + 100000*(50+noaa_number(SID));
nav[2]=SS + 100*MM + 10000*HH;
nav[3]=1;
{
int i;
int mm,dd;
int v=0;
mm=dd=0;
for(i=1;i<13;i++)
{
v+=days[i];
if(ED<=v)
{
mm=i;
dd=ED-v+days[i];
break;
}
}
nav[4]=dd+100*mm+10000*EY;
}
nav[5]=ES + 100*EM +10000*EH;
nav[6]=SMA/1000;
nav[7]=ECC/100;
nav[8]=INC/100;
nav[9]=MAN/100;
nav[10]=PER/100;
nav[11]=RAS/100;
nav[12]=2048;
nav[13]= 54128; /* FOV step angle */
nav[14]=ET%1000;
nav[45]=( south ? -1 : 1);
nav[46]=1;
nav[47]=MS;
nav[48]=167; /* scan period in ms */
nav[49]=0;
nav[50]=lines;
if(GAC) nav[50]*=3;
nav[51]=2048;
nav[52]=166667;
nav[53]=2500; /* step time */
}
#define I4(x) (i4(&head[x]))
static double C1 = 1.191044e-5;
static double C2 = 1.438769;
double c3[3];
double c2[6];
double c1[6];
double c0[6];
int intersect[3];
double p1[6];
double p2[6];
double p3[6];
double p4[6];
void load_cal_hdr(void)
{
unsigned char *p=head+RADCOEF;
double w;
int i;
p1[0] = ((double)i4(p))/1.e1;
p2[0] = ((double)i4(p+4))/1.e3;
p1[1] = ((double)i4(p+8))/1.e1;
p2[1] = ((double)i4(p+12))/1.e3;
p1[2] = ((double)i4(p+16))/1.e1;
p2[2] = ((double)i4(p+20))/1.e3;
w = i4(p+24);
w = w/1.e2;
p1[3] = C1*w*w*w;
p2[3] = C2*w;
p3[3] = i4(p+28);
p3[3] = p3[3]/1.e5;
p4[3] = i4(p+32);
p4[3] = p4[3]/1.e6;
w = i4(p+36);
w = w/1.e3;
p1[4] = C1*w*w*w;
p2[4] = C2*w;
p3[4] = i4(p+40);
p3[4] = p3[4]/1.e5;
p4[4] = i4(p+44);
p4[4] = p4[4]/1.e6;
w = i4(p+48);
w = w/1.e3;
p1[5] = C1*w*w*w;
p2[5] = C2*w;
p3[5] = i4(p+52);
p3[5] = p3[5]/1.e5;
p4[5] = i4(p+56);
p4[5] = p4[5]/1.e6;
#if 0
for(i=3;i<6;i++)
{
printf("%10.2f%10.4f%10f%10f\n",p1[i],p2[i],p3[i],p4[i]);
}
printf("\n");
#endif
sprintf((char *)cal,
"%10.2f%10.4f%10f%10f"
"%10.2f%10.4f%10f%10f"
"%10.2f%10.4f%10f%10f ",
p1[3],p2[3],p3[3],p4[3],
p1[4],p2[4],p3[4],p4[4],
p1[5],p2[5],p3[5],p4[5]);
}
void load_cal(unsigned char *x,char *output)
{
int i;
static int swit=0;
{
c0[0] = i4(x+CAL1); /* slope */
c0[0] = c0[0]/1.e7;
c1[0] = i4(x+CAL1+4); /* intercept */
c1[0] = c1[0]/1.e6;
c2[0] = i4(x+CAL1+8);
c2[0] = c2[0]/1.e7;
c3[0] = i4(x+CAL1+12);
c3[0] = c3[0]/1.e6;
intersect[0] = (c1[0] - c3[0])/(c2[0] - c0[0]);
if(swit) printf("%10.5f%10.5f%10.5f%10.5f %d\n",
c0[0],c1[0],c2[0],c3[0],intersect[0]);
c0[1] = i4(x+CAL2);
c0[1] = c0[1]/1.e7;
c1[1] = i4(x+CAL2+4);
c1[1] = c1[1]/1.e6;
c2[1] = i4(x+CAL2+8);
c2[1] = c2[1]/1.e7;
c3[1] = i4(x+CAL2+12);
c3[1] = c3[1]/1.e6;
intersect[1] = (c1[1] - c3[1])/(c2[1] - c0[1]);
if(swit) printf("%10.5f%10.5f%10.5f%10.5f %d\n",
c0[1],c1[1],c2[1],c3[1],intersect[1]);
c0[2] = i4(x+CAL3a);
c0[2] = c0[2]/1.e7;
c1[2] = i4(x+CAL3a+4);
c1[2] = c1[2]/1.e6;
c2[2] = i4(x+CAL3a+8);
c2[2] = c2[2]/1.e7;
c3[2] = i4(x+CAL3a+12);
c3[2] = c3[2]/1.e6;
intersect[2] = (c1[2] - c3[2])/(c2[2] - c0[2]);
if(i4(x+CAL3a+4)==0) intersect[2]=0;
if(swit) printf("%10.5f%10.5f%10.5f%10.5f %d\n",
c0[2],c1[2],c2[2],c3[2],intersect[2]);
c0[3]=i4(x+CAL3b);
c0[3]=c0[3]/1.e6;
c1[3]=i4(x+CAL3b+4);
c1[3]=c1[3]/1.e6;
c2[3]=i4(x+CAL3b+8);
c2[3]=c2[3]/1.e6;
if(swit) printf("%10.5f%10.5f%10.7f\n",c0[3],c1[3],c2[3]);
c0[4]=i4(x+CAL4);
c0[4]=c0[4]/1.e6;
c1[4]=i4(x+CAL4+4);
c1[4]=c1[4]/1.e6;
c2[4]=i4(x+CAL4+8);
c2[4]=c2[4]/1.e6;
if(swit) printf("%10.5f%10.5f%10.7f\n",c0[4],c1[4],c2[4]);
c0[5]=i4(x+CAL5);
c0[5]=c0[5]/1.e6;
c1[5]=i4(x+CAL5+4);
c1[5]=c1[5]/1.e6;
c2[5]=i4(x+CAL5+8);
c2[5]=c2[5]/1.e6;
if(swit) printf("%10.5f%10.5f%10.7f\n",c0[5],c1[5],c2[5]);
sprintf(output,
"%10.5f%10.5f%10.5f%10.5f"
"%10.5f%10.5f%10.5f%10.5f"
"%10.5f%10.5f%10.5f%10.5f"
"%10.5f%10.5f%10.7f"
"%10.5f%10.5f%10.7f"
"%10.5f%10.5f%10.7f ",
c0[0],c1[0],c2[0],c3[0],
c0[1],c1[1],c2[1],c3[1],
c0[2],c1[2],c2[2],c3[2],
c0[3],c1[3],c2[3],
c0[4],c1[4],c2[4],
c0[5],c1[5],c2[5]);
}
swit=0;
}
/*
* construct a machine-ordered integer out of 4 big-endian bytes
*/
int i4(unsigned char *x)
{
int j;
if(x[0]&0x80)
{
j=x[3]^0xff;
j=j+256*(x[2]^0xff);
j=j+256*256*(x[1]^0xff);
j=j+256*256*256*((x[0]&0x7f)^0x7f);
j=-j;
}
else
{
j=x[3]+256*(x[2]+256*(x[1]+256*x[0]));
}
return j;
}
/*
* return noaa number from
* integer stored in file
*/
int noaa_number(int i)
{
if(i==6) return 8;
if(i==7) return 9;
if(i==8) return 10;
if(i==1) return 11; /* warning: before 87, this was noaa5 */
if(i==5) return 12;
/* 13 didn't make it. dead on arrival */
if(i==3) return 14;
if(i==4) return 15; /* was noaa7 earlier */
return -1;
}
------------------------------ conlev1b.c ------------------------------------
Tom
--
+-----------------------------------------------------------------------------+
* Tom Yoksas UCAR Unidata Program *
* (303) 497-8642 (last resort) P.O. Box 3000 *
* address@hidden Boulder, CO 80307 *
* Unidata WWW Service http://www.unidata.ucar.edu/*
+-----------------------------------------------------------------------------+