#Original Function Location: /themis/lib/dav_lib/public_library/latest/crism_science.dvrc
/themis/lib/dav_lib/library/crism_science.dvrc
Source Code for Function: "crism_destripe()" in "crism_science.dvrc" (Public)
crism_science_version=1.12
#process_crism - modified by c.edwards, adapted from l.mayorga 3/09**
#get_crism - modified by c.edwards, adapted from j.bandfield 3/09 taken from crism.dav*
#crism_atm - modified by c.edwards, adapted from j.bandfield 3/09 taken from crism.dav*
#crism_resample - modified by c.edwards, adapted from j.bandfield 3/09 taken from crism.dav
#crism_destripe - created by c.edwards, modified by j.bandfield 3/09 taken from crism.dav
#project_crism - created by c.edwards, modified by j.bandfield 3/09 taken from crism.dav
#project_crism_index - modified by c.edwards, created by l.mayorga 3/09 taken from crism.dav**
#create_indices - modified by c.edwards, created by j.bandfield 3/09 taken from crism.dav
#create_indices_key - modified by c.edwards, created by j.bandfield 3/09 taken from crism.dav
#crism_speclib - modified by c.edwards, created by j.bandfield 3/09 taken from crism.dav
#get_crism_list - modified by c.edwards, created by j.bandfield 3/09 taken from crism.dav
#get_proj_info - created by j.bandfield 2/10
#* modified by l.mayorga 7/16/09
#** modified by l.mayorga 7/23/09
# modified to use the appserver to get crism data.
# "fixed" destripe/atm correction 10/7/09 - j.bandfield/c.edwards
#
# added option to process_crism() to allow the user to not roate the data to ~N-up
# this allows for comparisons to other crism publications - c.edawrds 11-13-10
#
# modified get_crism() to use the pds source file by default...this is a new option in get_image (misc.dvrc)
# c.edwards 8-25-2011
#
# modified get_crism() to use the new pds url
# c.edwards 9-14-2011
# modified project_crism() to project to a set resolution for mosaicking added moon option
# modified get_proj_info() to better parse gdalinfo output
# modified get_crism() and process_crism() to find specific atm files for atm. correction
# now handles trr3 images other minor fixes
# e. amador, c. edwards, j. bandfield 9-7-11
define process_crism(raw, resample, atm, sl, index, destripe,wget,rotate,src) {
#tested cse 3/30/09
if($ARGC!=1) {
printf("Download and process any CRISM images, short and long, radiance or I/F\n\n")
printf("$1=product id (truncated from full ID)\n\n")
printf("Optional:\n")
printf("\tsrc=download from the pds source and not the MSFF (Default=0)\n")
printf("\tsl=process both the short and long data product (Default=0, long only)\n")
printf("\traw=return just the output structure (other options will not work with this, Default=0)\n")
printf("\tresample=resample the image to a common scale, uses crism_resample() (the detector is not wavelength uniform, Default=0)\n")
printf("\tatm=do the atmospheric correction for the long products, uses crism_atm() (Default=1)\n")
printf("\tdestripe=destripe and destripe the data, uses crism_destripe() (Default=1)\n")
printf("\tindex=generate standard indices + a key as a structure element (Default=1)\n")
printf("\trotate=rotate the image to ~North up (Default=1)\n")
printf("modified by c.edwards, adapted from l.mayorga, 3/9\n")
return(null)
}
product=$1
if (HasValue(raw)==0) raw=0
if (HasValue(resample)==0) resample=0
if (HasValue(atm)==0) atm=1
if (HasValue(sl)==0) sl=0
if (HasValue(index)==0) index=1
if (HasValue(destripe)==0) destripe=1
if (HasValue(rotate)==0) rotate=1
if (HasValue(url)==0) url=0
if (HasValue(src)==0) src=0
if(HasValue(wget)==0) wget=0
if(atof(version()[18:])>=2.0) {
copy=0
} else {
copy=1
}
if (wget==1) {
copy=1
}
#get the data
data=get_crism(product,raw=raw,sl=sl,wget=copy,src=src)
if(raw==1) {
return(data)
}
# Check to see if this is a TRR2 or TRR3 image, if TRR3 destripe=0
tmpdir=sprintf("%s/%s",$DV_CRISM_CACHE,product)
if (HasValue(syscall("ls "+tmpdir+"/*_TRR3.IMG")) ==1) {
destripe=0
}
if(rotate==1) {
#flip the data to orient it with ~N up
data.long.data=translate(translate(data.long.data,from=x,to=x,flip=1),y,y,flip=1)
data.long.anc.ddr=translate(translate(data.long.anc.ddr,from=x,to=x,flip=1),y,y,flip=1)
data.long.wl_array=translate(translate(data.long.wl_array,from=x,to=x,flip=1),y,y,flip=1)
data.long.at_array=translate(translate(data.long.at_array,from=x,to=x,flip=1),y,y,flip=1)
}
#set the null value
data.long.data [where data.long.data==65535]=0
if(sl==1) {
data.short.data [where data.short.data==65535]=0
}
verbose=0
if(atm==1) data=crism_atm(data=data)
if(resample==1) data=crism_resample(data=data)
if(destripe==1) data.long.data=crism_destripe(data=data.long.data,v=1)
if(destripe==2) data.long.data=crism_destripe_old(data=data.long.data)
if(index==1) {
data.long.index=create_indices(data=data.long.data)
data.long.index[where (data.long.data[,,50]==0)]=0
data.long.index_key=create_indices_key(data.long.index)
if(sl==1) {
data.short.index=create_indices(data=data.short.data)
data.long.index[where (data.short.data[,,10]==0)]=0
data.short.index_key=create_indices_key(data.short.index)
}
}
verbose=3
return(data)
}
define get_crism(ra,raw,sl,force,verb,ignore,set_ignore,wget,src) {
#tested-cse 3/30/09
if($ARGC!=1) {
printf("Download any CRISM image, short and long, radiance or I/F\n\n")
printf("$1=product id (truncated from full ID)\n\n")
printf("Optional:\n")
printf("\tsrc=download from the pds source and not ASU (Default=1)\n")
printf("\tsl=process both the short and long data product (Default=0, long only)\n")
printf("\traw=return just the output structure (other options will not work with this, Default=0)\n")
printf("\tforce=force the redownload of products (Default=0), uses $DV_CRISM_CACHE\n")
printf("\tverb=set the verbosity (Default=0)\n")
printf("\tignore=the ignore value (Default=0)\n\n")
printf("\tset_ignore=set the value of value of the image (Default=1)\n")
printf("\tsrc=download the data from the PDS rather than ASU (default=1)\n")
printf("\twget = force the use of wget (Default=0)\n\n")
printf("modified by c.edwards, adapted from j.bandfield 3/09\n")
return(null)
}
product=$1
global(verbose)
verbold=verbose
if(HasValue(src)==0) src=1
if(HasValue(wget)==0) wget=0
if(atof(version()[18:])>=2.0) {
copy=1
} else {
copy=0
}
if (wget==1) {
copy=0
}
printf("Copy %i\n",copy)
printf("Running get CRISM...")
if (HasValue(raw)==0) raw=0
if (HasValue(sl)==0) sl=0
if (HasValue(verb)==0) verbose=0
if (HasValue(force)==0) force=0
if (HasValue(ignore)==0) ignore=0
if (HasValue(set_ignore)==0) set_ignore=0
tmpdir=sprintf("%s/%s",$DV_CRISM_CACHE,product)
printf("%s\n",tmpdir)
system(sprintf("mkdir -p %s",tmpdir))
system(sprintf("chmod 777 %s",tmpdir))
exists=syscall("ls "+tmpdir)
if(HasValue(exists)) {
Strue=length(grep(exists,"S_TRR"))==2
Ltrue=length(grep(exists,"S_TRR"))==2
SDtrue=length(grep(exists,"S_DDR"))==2
LDtrue=length(grep(exists,"S_DDR"))==2
} else {
Strue=0
Ltrue=0
SDtrue=0
LDtrue=0
}
if(sl!=0) {
if(Strue==1) {
url_trrs_lbl=syscall("ls "+tmpdir+"/*S_TRR?.LBL")[,1]
url_trrs_img=syscall("ls "+tmpdir+"/*S_TRR?.IMG")[,1]
} else {
url_trrs_lbl=get_image(product+"_S",instrument="crism",type="TRDR-LBL",src=src)
url_trrs_img=get_image(product+"_S",instrument="crism",type="TRDR-IMG",src=src)
}
if(copy==1) {
trrs_filename=basename(load_pds(url_trrs_lbl,data=0).product_id[2:],"\"")
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_trrs_lbl,$TMPDIR+"/trrs.lbl"))
trrs_filename=basename(load_pds($TMPDIR+"/trrs.lbl",data=0).product_id[2:],"\"")
}
if(trrs_filename=="\"") {
printf("Could not download the TRR %s_S\n",product)
return(null)
}
trrs_filename=tmpdir+"/"+trrs_filename
if(fexists(trrs_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",trrs_filename)
if(copy==1) {
copy(url_trrs_lbl,trrs_filename+".LBL",force=force)
copy(url_trrs_img,trrs_filename+".IMG",force=force)
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_trrs_lbl,trrs_filename+".LBL"))
system(sprintf("%s '%s' -O %s","wget -nd",url_trrs_img,trrs_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",trrs_filename)
}
if (SDtrue==1) {
url_ddrs_lbl=syscall("ls "+tmpdir+"/*S_DDR?.LBL")[,1]
url_ddrs_img=syscall("ls "+tmpdir+"/*S_DDR?.IMG")[,1]
} else {
url_ddrs_lbl=get_image(product+"_S",instrument="crism",type="DDR-LBL",src=src)
url_ddrs_img=get_image(product+"_S",instrument="crism",type="DDR-IMG",src=src)
}
if(copy==1) {
ddrs_filename=basename(load_pds(url_ddrs_lbl,data=0).product_id[2:],"\"")
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_ddrs_lbl,$TMPDIR+"/ddrs.lbl"))
ddrs_filename=basename(load_pds($TMPDIR+"/ddrs.lbl",data=0).product_id[2:],"\"")
}
if(ddrs_filename=="\"") {
printf("Could not download the DDR %s_S\n",product)
return(null)
}
ddrs_filename=tmpdir+"/"+ddrs_filename
if(fexists(ddrs_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",ddrs_filename)
if(copy==1) {
copy(url_ddrs_lbl,ddrs_filename+".LBL",force=force)
copy(url_ddrs_img,ddrs_filename+".IMG",force=force)
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_ddrs_lbl,ddrs_filename+".LBL"))
system(sprintf("%s '%s' -O %s","wget -nd",url_ddrs_img,ddrs_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",ddrs_filename)
}
}
if (Ltrue==1) {
url_trrl_lbl=syscall("ls "+tmpdir+"/*L_TRR?.LBL")[,1]
url_trrl_img=syscall("ls "+tmpdir+"/*L_TRR?.IMG")[,1]
} else {
url_trrl_lbl=get_image(product+"_L",instrument="crism",type="TRDR-LBL",src=src)
url_trrl_img=get_image(product+"_L",instrument="crism",type="TRDR-IMG",src=src)
}
if(copy==1) {
trrl_filename=basename(load_pds(url_trrl_lbl,data=0).product_id[2:],"\"")
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_trrl_lbl,$TMPDIR+"/trrl.lbl"))
trrl_filename=basename(load_pds($TMPDIR+"/trrl.lbl",data=0).product_id[2:],"\"")
}
if(trrl_filename=="\"") {
printf("Could not download the TRR %s_L\n",product)
return(null)
}
trrl_filename=tmpdir+"/"+trrl_filename
if(fexists(trrl_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",trrl_filename)
if(copy==1) {
copy(url_trrl_lbl,trrl_filename+".LBL",force=force)
copy(url_trrl_img,trrl_filename+".IMG",force=force)
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_trrl_lbl,trrl_filename+".LBL"))
system(sprintf("%s '%s' -O %s","wget -nd",url_trrl_img,trrl_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",trrl_filename)
}
if (LDtrue ==1) {
url_ddrl_lbl=syscall("ls "+tmpdir+"/*L_DDR?.LBL")[,1]
url_ddrl_img=syscall("ls "+tmpdir+"/*L_DDR?.IMG")[,1]
} else {
url_ddrl_lbl=get_image(product+"_L",instrument="crism",type="DDR-LBL",src=src)
url_ddrl_img=get_image(product+"_L",instrument="crism",type="DDR-IMG",src=src)
}
if(copy==1) {
ddrl_filename=basename(load_pds(url_ddrl_lbl,data=0).product_id[2:],"\"")
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_ddrl_lbl,$TMPDIR+"/ddrl.lbl"))
ddrl_filename=basename(load_pds($TMPDIR+"/ddrl.lbl",data=0).product_id[2:],"\"")
}
if(ddrl_filename=="\"") {
printf("Could not download the DDR %s_L\n",product)
return(null)
}
ddrl_filename=tmpdir+"/"+ddrl_filename
if(fexists(ddrl_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",ddrl_filename)
if(copy==1) {
copy(url_ddrl_lbl,ddrl_filename+".LBL",force=0)
copy(url_ddrl_img,ddrl_filename+".IMG",force=0)
} else {
system(sprintf("%s '%s' -O %s","wget -nd",url_ddrl_lbl,ddrl_filename+".LBL"))
system(sprintf("%s '%s' -O %s","wget -nd",url_ddrl_img,ddrl_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",ddrl_filename)
}
data=struct(long,ddrl)
data.ddrl=load_pds(ddrl_filename+".LBL")
data.long=load_pds(trrl_filename+".LBL",data=0)
data.long.file.image=read(trrl_filename+".LBL")
if (sl!=0) {
add_struct(data,name="short")
add_struct(data,name="ddrs")
data.ddrs=load_pds(ddrs_filename+".LBL")
data.short=load_pds(trrs_filename+".LBL",data=0)
data.short.file.image=read(trrs_filename+".LBL")
}
# Get wavelength arrays
if(HasValue(get_struct(data.long,"mro\:wavelength_file_name"))) {
wll_filename=basename(get_struct(data.long,"mro\:wavelength_file_name")[2:],".IMG\"")
} else {
wll_filename=basename(get_struct(data.long,"mro_wavelength_file_name")[2:],".IMG\"")
}
# surprise they changed the path again
wll_lbl="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wll_filename//".LBL"
wll_img="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wll_filename//".IMG"
# wll_lbl="http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wll_filename//".LBL"
# wll_img="http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wll_filename//".IMG"
if(fexists(tmpdir+"/"+wll_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",tmpdir+"/"+wll_filename)
if(copy==1) {
copy(wll_lbl,tmpdir+"/"+wll_filename+".LBL",force=force)
copy(wll_img,tmpdir+"/"+wll_filename+".IMG",force=force)
} else {
system(sprintf("%s %s -O %s","wget -nd", wll_lbl,tmpdir+"/"+wll_filename+".LBL"))
system(sprintf("%s %s -O %s","wget -nd", wll_img,tmpdir+"/"+wll_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",tmpdir+"/"+wll_filename)
}
data.long.wl_array=load_pds(tmpdir+"/"+wll_filename+".LBL")
if (sl!=0) {
if(HasValue(get_struct(data.long,"mro\:wavelength_file_name"))) {
wls_filename=basename(get_struct(data.short,"mro\:wavelength_file_name")[2:],".IMG\"")
} else {
wls_filename=basename(get_struct(data.short,"mro_wavelength_file_name")[2:],".IMG\"")
}
# surprise they changed the path again
wls_lbl="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wls_filename//".LBL"
wls_img="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wls_filename//".IMG"
# wls_lbl="http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wls_filename//".LBL"
# wls_img="http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wls_filename//".IMG"
if(fexists(tmpdir+"/"+wls_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",tmpdir+"/"+wls_filename)
if(copy==1) {
copy(wls_lbl,tmpdir+"/"+wls_filename+".LBL",force=force)
copy(wls_img,tmpdir+"/"+wls_filename+".IMG",force=force)
} else {
system(sprintf("%s %s -O %s","wget -nd",wls_lbl,tmpdir+"/"+wls_filename+".LBL"))
system(sprintf("%s %s -O %s","wget -nd",wls_img,tmpdir+"/"+wls_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",tmpdir+"/"+wls_filename)
}
data.short.wl_array=load_pds(tmpdir+"/"+wls_filename+".LBL")
}
# Get atmospheric transmission arrays
# Use these for the atmospheric correction (long wave only)
# find the text string that contains the spacecraft clock start count and converts it from text to double format, sets to the varibale sc_time
sc_time=atod(data.long.spacecraft_clock_start_count[4:13])
#Get second part of AT filename from the Wavelength Filename in the TRR2 file.
if(HasValue(get_struct(data.long,"mro\:wavelength_file_name"))) {
wl_filename3=get_struct(data.long,"mro\:wavelength_file_name")
} else {
wl_filename3=get_struct(data.long,"mro_wavelength_file_name")
}
#the varibale at_end is a text string of the last 7 digits to be used in the AT_Filename ex) CDR420843667218_AT"-->0000000L<--"_7.LBL
at_end=wl_filename3[20:26]
#Getting CDR Index Tab from the PDS website
#download the CDR Index Tab from the PDS website and set to variable cdrindex
cdrurl="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/index/cdrstatic_index.tab"
if(copy==1) {
copy(cdrurl,$TMPDIR+"/cdrstatic_index.tab",force=1)
} else {
cdrindex="wget -nd "+cdrurl
# create the string "cd $TMPDIR; cdrindex
command=sprintf("cd %s; %s",$TMPDIR,cdrindex)
#open the directory $TMPDIR and run the run the variable cdrindex which wget to download the CDR Indext temporarily into the $TMPDIR
system(command)
}
#look through the cdrstatic_index.tab located in $TMPDIR and grap all lines that contain "CDR/AT/" store as a text files into the variable atfiles
atfiles=syscall("cat "+$TMPDIR+"/cdrstatic_index.tab | grep \"CDR/AT/\"")
#Converts text to double in atfiles, characters 22 through 31, this is the start sclk time for the CDR files, new file called atfiles2
atfiles2=atod(atfiles[22:31])
# Loop through all the sclk times and create a file with each sclk time removing duplicates and saving to file outatfiles2, outatfiles saves the same lines to be used later when putting together full filename
outatfiles2=atfiles2[,1]
outatfiles=atfiles[,1]
for (i=1; i<=(dim(atfiles2)[2])-1; i+=1) {
if (atfiles2[,i] != atfiles2[,(i+1)]) {
outatfiles2=cat(outatfiles2,atfiles2[,(i+1)],axis=y)
outatfiles=cat(outatfiles,atfiles[,(i+1)],axis=y)
}
}
# Comparing sc_time to the list of available CDR files
if (outatfiles2[,1]>sc_time) {
outsclk=outatfiles2[,1]
line=outatfiles[,1]
}
numtimes=dim(outatfiles2)[2]
if (outatfiles2[,numtimes]
outsclk=outatfiles2[,numtimes]
line=outatfiles[,numtimes]
}
for (i=1; i<=numtimes-1; i+=1) {
if (outatfiles2[,i]<=sc_time && outatfiles2[,i+1]>sc_time) {
outsclk=outatfiles2[,i]
line=outatfiles[,i]
}
}
at_filename1=line[17:34]
version_number=line[44]
#Atmsospheric transmission array filename
atl_filename=sprintf("%s%s%s%s",at_filename1,at_end,"L_",version_number)
#atl_filename=basename(get_struct(data.long,"mro\:wavelength_file_name")[2:],".IMG\"")
#atl_filename=atl_filename[1:5]//"0000000000_AT"//wll_filename[19:]
# moved file location...surprise!!!
atl_lbl="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//atl_filename//".LBL"
atl_img="http://geodata.rsl.wustl.edu/mro/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//atl_filename//".IMG"
# atl_lbl="http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//atl_filename//".LBL"
# atl_img="http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//atl_filename//".IMG"
if(fexists(tmpdir+"/"+atl_filename+".LBL")==0 || force) {
printf("Downloading: %s\n",tmpdir+"/"+atl_filename)
if(copy==1) {
copy(atl_lbl,tmpdir+"/"+atl_filename+".LBL",force=force)
copy(atl_img,tmpdir+"/"+atl_filename+".IMG",force=force)
} else {
system(sprintf("%s %s -O %s","wget -nd",atl_lbl,tmpdir+"/"+atl_filename+".LBL"))
system(sprintf("%s %s -O %s","wget -nd",atl_img,tmpdir+"/"+atl_filename+".IMG"))
}
} else {
printf("File Exists: %s\n",tmpdir+"/"+atl_filename)
}
data.long.at_array=load_pds(tmpdir+"/"+atl_filename+".LBL")
verbose=verbold
if (raw==1) {
return(data)
}
if (raw!=1) {
out=struct(long)
out.long=struct(anc,wl_array,at_array,id,data)
out.long.anc=struct(ls,solar_dist,ddr,ddrkey)
out.long.anc.summing=data.long.pixel_averaging_width
out.long.data=remove_struct(data.long.file,"image")
if(set_ignore==1) {
out.long.data[where out.long.data==65535]=ignore
}
out.long.id=data.long.product_id[,1]
out.long.wl_array=remove_struct(data.long.wl_array.file.image,"data")
out.long.at_array=remove_struct(data.long.at_array.file.image,"data")
out.long.anc.ls=data.ddrl.solar_longitude
out.long.anc.solar_dist=data.ddrl.solar_distance
out.long.anc.ddrkey=data.ddrl.file.image.band_name
out.long.anc.ddr=remove_struct(data.ddrl.file.image,"data")
if (sl!=0) {
out.short=struct(anc,wl_array,id,data)
out.short.anc=struct(ls,solar_dist,ddr,ddrkey)
out.short.anc.summing=data.short.pixel_averaging_width
out.short.data=remove_struct(data.short.file,"image")
if(set_ignore==1) {
out.short.data[where out.short.data==65535]=ignore
}
out.short.id=data.short.product_id[,1]
out.short.wl_array=remove_struct(data.short.wl_array.file.image,"data")
out.short.anc.ls=data.ddrs.solar_longitude
out.short.anc.solar_dist=data.ddrs.solar_distance
out.short.anc.ddrkey=data.ddrs.file.image.band_name
out.short.anc.ddr=remove_struct(data.ddrs.file.image,"data")
}
}
printf("Done.\n")
return(out)
}
define crism_atm(data) {
#tested-cse 3/30/09
if(HasValue(data)==0 && $ARGC==1) {
data=$1
} else if (HasValue(data)==0 && $ARGC==0) {
printf("Perform the atmospheric correction for long products\n")
printf("$1=standard CRISM structure (long or long/short) produced by get_crism()\n")
printf("Or\n")
printf("data=standard CRISM structure (long or long/short) produced by get_crism()\n")
printf("modified by c.edwards, created by j.bandfield, 3/9\n")
return(null)
}
printf("Running CRISM atmospheric correction...")
data.long.at_array=convolve2(data.long.at_array,clone(1.,50,1,1),ignore=65535.0)
#Modifed by l.mayorga 7/16/09
summing = data.long.anc.summing
if (dim(data.long.data)[3]==438){
R1980 = avg(data.long.data[,,288:291],axis=z)
R2007 = avg(data.long.data[,,284:287],axis=z)
A1980 = avg(data.long.at_array[,,288:291],axis=z)
A2007 = avg(data.long.at_array[,,284:287],axis=z)
duplicate = 437
} else {
R1980 = avg(data.long.data[,,30],axis=z)
R2007 = avg(data.long.data[,,29],axis=z)
A1980 = avg(data.long.at_array[,,30],axis=z)
A2007 = avg(data.long.at_array[,,29],axis=z)
duplicate = 54
}
expon = ln(R1980/R2007)/ln(A1980/A2007)
if (summing==1) {
for (i=8; i<=611; i+=1) {
data.long.data[i,,2:] = data.long.data[i,,2:] / data.long.at_array[i,,2:]^clone(expon[i],z=duplicate)
}
} else if (summing==2) {
for (i=5; i<=305; i+=1) {
data.long.data[i,,2:] = data.long.data[i,,2:] / data.long.at_array[i,,2:]^clone(expon[i],z=duplicate)
}
} else if (summing==5 || summing==10) {
for (i=(10/summing+1); i<=(640/summing)-(30/summing); i+=1) {
data.long.data[i,,2:] = data.long.data[i,,2:] / data.long.at_array[i,,2:]^clone(expon[i],z=duplicate)
}
}
printf("Done.\n")
}
define crism_resample(data) {
#tested cse 3/30/09 --needs more
if(HasValue(data)==0 && $ARGC==1) {
data=$1
} else if (HasValue(data)==0 && $ARGC==0) {
printf("Resample the image to be uniform across the entire detector\n")
printf("$1=standard CRISM structure (long or long/short) produced by get_crism()\n")
printf("Or\n")
printf("data=standard CRISM structure (long or long/short) produced by get_crism()\n\n")
printf("modified by c.edwards, created by j.bandfield, 3/9\n")
return(null)
}
printf("Running CRISM resample...")
ave_wll=avg(data.long.wl_array[int(260/data.long.anc.summing):int(359/data.long.anc.summing)],axis=x)
for (i=(30/data.long.anc.summing); i<=(634/data.long.anc.summing-1); i+=1) {
data.long.data[i,,2:]=translate(resample(translate(data.long.data[i,,2:],z,z,flip=1),translate(data.long.wl_array[i,,2:],z,z,flip=1),translate(ave_wll[,,2:],z,z,flip=1)).data,z,z,flip=1)
}
data.long.wl_array=ave_wll
verbose=0
if (HasValue(data.short)!=0) {
ave_wls=avg(data.short.wl_array[int(270/data.long.anc.summing):int(369/data.long.anc.summing)],axis=x)
for (i=(26/data.long.anc.summing); i<=(626/data.long.anc.summing-1); i+=1) {
data.short.data[i]=translate(resample(translate(data.short.data[i],z,z,flip=1),translate(data.short.wl_array[i],z,z,flip=1),translate(ave_wls,z,z,flip=1)).data,z,z,flip=1)
}
data.short.wl_array=ave_wls
}
verbose=1
printf("Done.\n")
return(data)
}
define crism_destripe(r_size,f_size,ignore,sigmamult,data,v) {
#tested cse 3/30/09
if(HasValue(r_size)==0) r_size=21
if(HasValue(z_size)==0) z_size=9
if(HasValue(f_size)==0) f_size=3
if(HasValue(ignore)==0) ignore=0
if(HasValue(sigmamult)==0) sigmamult=1.0
if(HasValue(v)==0) v=0
if(HasValue(data)==0 && $ARGC==1) {
data=$1
} else if (HasValue(data)==0 && $ARGC==0) {
printf("Perform the despike/destripe for either long or short products\n")
printf("$1=data element from the standard crism structure (long or short) produced by get_crism()\n")
printf("Or\n")
printf("data=data element from the standard crism structure (long or short) produced by get_crism()\n\n")
printf("Options:\n")
printf("\tr_size=residual column noise filter size (Default=21)\n")
printf("\tf_size=spike detection filter size (Default=3)\n")
printf("\tz_size=spectral filter size for replacement values (Default=9)\n")
printf("\tsigmamult=number of standard deviations from the mean for spike detection (Default=1.0)\n")
printf("\tignore=null value for the image (Default=0)\n\n")
printf("modified by j.bandfield/c.edwards, 10/09\n")
return(null)
}
printf("Running CRISM destripe...")
#construct the mask
mask=byte(data[,,1]*0)
mask[where (data[,,5]!=ignore)]=1
for(j=1;j<=dim(data)[2];j+=1) {
if(v>0) printf("Line: %i\n",j)
#extract the line and replacement array
line=data[,j,5:]
#create the smoothed array and calculate the average+stddev
con=convolve2(line,clone(1.,f_size,1,1),ignore=ignore)
avg=avg(line-con,both=1,ignore=ignore,axis=x)
#deal with the low> spike and deal with the high spikes
line[where ((-(avg.stddev)*sigmamult)>(line-con))]=ignore
line[where (((avg.stddev)*sigmamult)<(line-con))]=ignore
sigma=translate(boxfilter(translate(line,z,y),x=1,y=z_size,ignore=ignore,verbose=1).sigma,y,z)
#filter the data in spectral direction
line2=translate(boxfilter(translate(line,z,y),x=1,y=z_size,ignore=ignore),y,z)
#replace the original data where spikes were detected
line[where line==ignore]=line2
#replace the original data where no good replacement data could be created
line[where line==ignore]=data[,j,5:]
line[where sigma==0]=data[,j,5:]
#back in to the original array for memory's sake
data[,j,5:]=float(line)
}
#reset the ignore values
data[where mask==0]=ignore
#run deplaid last (seems to make the noise a bit clearer)
#this option usees a directional deplaid for older versions of davinci
if(dim(data)[3]==438 && atof(version()[18:])>=2.0) {
for (i=5; i<=168; i+=20) {
if (i<149) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0,axis=y)
} else {
data[,,i:168]=thm.deplaid(data[,,i:168],ignore=ignore,b10=0,axis=y)
}
}
for (i=187; i<=438; i+=20) {
if (i<419) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0,axis=y)
} else {
data[,,i:438]=thm.deplaid(data[,,i:438],ignore=ignore,b10=0,axis=y)
}
}
}
#this option uses the standard deplaid if the directional version is not available
if(dim(data)[3]==438 && atof(version()[18:])<2.0) {
for (i=5; i<=168; i+=20) {
if (i<149) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0)
} else {
data[,,i:168]=thm.deplaid(data[,,i:168],ignore=ignore,b10=0)
}
}
for (i=187; i<=438; i+=20) {
if (i<419) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0)
} else {
data[,,i:438]=thm.deplaid(data[,,i:438],ignore=ignore,b10=0)
}
}
}
#remove residual line noise
avg=avg(data,ignore=ignore,axis=y)
con=convolve2(avg,clone(1.,r_size,1,1),ignore=ignore)
data=data-float((avg-con))
printf("done.\n")
return(data)
}
define crism_destripe_old(r_size,f_size,ignore,sigmamult,data,v) {
#tested cse 3/30/09
if(HasValue(r_size)==0) r_size=7
if(HasValue(z_size)==0) z_size=9
if(HasValue(f_size)==0) f_size=3
if(HasValue(ignore)==0) ignore=0
if(HasValue(sigmamult)==0) sigmamult=1.25
if(HasValue(v)==0) v=0
if(HasValue(data)==0 && $ARGC==1) {
data=$1
} else if (HasValue(data)==0 && $ARGC==0) {
printf("Perform the despike/destripe for either long or short products\n")
printf("$1=data element from the standard crism structure (long or short) produced by get_crism()\n")
printf("Or\n")
printf("data=data element from the standard crism structure (long or short) produced by get_crism()\n\n")
printf("Options:\n")
printf("\tr_size=replacement image filter size (Default=5)\n")
printf("\tf_size=spike detection filter size (Default=3)\n")
printf("\tsigmamult=number of standard deviations from the mean for spike detection (Default=1.25)\n\n")
printf("\tignore=null value for the image (Default=0)\n\n")
printf("modified by c.edwards, created by j.bandfield, 3/9\n")
return(null)
}
#construct the mask
mask=byte(data[,,1]*0)
mask[where (data[,,5]!=ignore)]=1
# plim(0,700,0.02,.25)
for(k=1;k<=dim(data)[3];k+=1) {
printf("%i\n",k)
#extract a single band to work with (deplaid takes care of band correlated noise)
#addend/beginning fixes
start=k-(int(z_size/2))
end=k+(int(z_size/2))
if(start<1) start=1
if(end>dim(data)[3]) end=dim(data)[3]
replace=data[,,start:end]
#get a windowed average and stddev
pxavg=boxfilter(replace,x=r_size,ignore=ignore,verbose=1)
replace[where mask==0]=ignore
#construct the replace array from the windowed average/stddev arrays
replace[where replace pxavg.mean+pxavg.sigma]=ignore
middle=k-start+1
replace[,,middle]=ignore
replace=float(avg(replace,axis=z,ignore=ignore))
#convlove the array to make a uniform scene
replace2=boxfilter(replace,size=3,ignore=ignore)
replace[where replace==ignore]=replace2
replace[where mask==0]=ignore
#clean up the replace array
avg=avg(replace,axis=y,ignore=ignore)
con=convolve2(avg,clone(1.,51,1,1),ignore=ignore)
replace=replace-float(avg-con)
for(j=1;j<=dim(data)[2];j+=1) {
#extract the line and replacement array
line_old=data[,j,k]
line=line_old
replace_line=replace[,j]
#create the smoothed array and calculate the average+stddev
#con_bf=boxfilter(line,size=f_size,ignore=ignore)
con=convolve2(line,clone(1.,f_size,f_size,1),ignore=ignore)
avg=avg(line-con,both=1,ignore=ignore)
#deal with the low> spike
line[where ((-(avg.stddev)*sigmamult)>(line-con))]=replace_line
#deal with the high spikes
line[where (((avg.stddev)*sigmamult)<(line-con))]=replace_line
data[,j,k]=float(line)
#xplot(line_old,line,ignore=0)
}
}
#reset the ignore values
data[where mask==0]=ignore
#run deplaid first (seems to make the noise a bit clearer)
# data=thm.deplaid(data,ignore=ignore,b10=0)
if(dim(data)[3]==438 && atof(version()[18:])>=2.0) {
for (i=4; i<=168; i+=20) {
if (i<149) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0,axis=y)
} else {
data[,,i:168]=thm.deplaid(data[,,i:168],ignore=ignore,b10=0,axis=y)
}
}
for (i=187; i<=438; i+=20) {
if (i<419) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0,axis=y)
} else {
data[,,i:438]=thm.deplaid(data[,,i:438],ignore=ignore,b10=0,axis=y)
}
}
}
if(dim(data)[3]==438 && atof(version()[18:])<2.0) {
for (i=4; i<=168; i+=20) {
if (i<149) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0)
} else {
data[,,i:168]=thm.deplaid(data[,,i:168],ignore=ignore,b10=0)
}
}
for (i=187; i<=438; i+=20) {
if (i<419) {
data[,,i:(i+20)]=thm.deplaid(data[,,i:(i+20)],ignore=ignore,b10=0)
} else {
data[,,i:438]=thm.deplaid(data[,,i:438],ignore=ignore,b10=0)
}
}
}
#remove residual line noise
avg=avg(data,ignore=ignore,axis=y)
con=convolve2(avg,clone(1.,r_size,1,1),ignore=ignore)
data=data-float((avg-con))
return(data)
}
define project_crism(proj,planet,ddr,clat,clon,outstruct,sampling,ignore,manual) {
#tested cse 03/30/09
#updated jlb 09/07/11
if($ARGC!=2) {
printf("\nProject a crism product using gdal from control points in the DDR image\n")
printf("$1 = the CRISM DDR data array, 2 band lat/lon array\n")
printf("$2 = fully proccessed CRISM image (one or three band byte)\n")
printf("Options:\n")
printf("\tproj=projection type (\"sinu\",\"simp\", Default=\"simp\"\n")
printf("\tplanet=1-use the Mars ellipsoid reference 2-Earth epsg:4326 3-Moon (Default=1)\n")
printf("\tddr=use the full ddr array (will choose the right bands, Default=1)\n")
printf("\tsampling=pixel sampling in meters (default is native sampling)\n")
printf("\tclon=use a user specified center longitude for \"sinu\" projection, Default is center of image\n")
printf("\tclat=use a user specified center latitude for \"simp\" projection, Default is 0\n")
printf("\toutstruct=run gdalinfo and return geometry info, Default is 1\n")
printf("\tmanual=provide list of control points instead fo lat/lon backplanes\n")
printf("\t(as $1 with 4 columns x,y,lon,lat), Default is 0\n")
printf("\tignore=data values treated as null by GDAL, default is 0\n")
printf("NOTE: requires gdalwarp and gdal_translate >1.5.1, specified by $DV_GDAL_PATH\n\n")
printf("modified by j.bandfield,c.edwards, written by c.edwards 3/09\n\n")
return(null)
}
if(HasValue(proj)==0) proj="sinu"
if(HasValue(planet)==0) planet=1
if(HasValue(ddr)==0) ddr=1
if(HasValue(outstruct)==0) outstruct=1
if(HasValue(sampling)==0) sampling=0
if(HasValue(ignore)==0) ignore=0
if(HasValue(manual)==0) manual=0
latlon=$1
img=$2
printf(""+ddr+"")
printf(""+planet+"")
if (sampling==0) {
sampopt=" "
} else {
sampopt="-tr "+sampling+" "+sampling
}
if(ddr==1 && manual==0) {
latlon=latlon[,,4:5]
}
dlatlon=dim(latlon)
nbands=dim(img)[3]
for (bnum=1; bnum<=nbands; bnum+=1) {
write(img[,,bnum],$TMPDIR+"/project_crism_temp_image.tif",type=tif,force=1)
if(fexists($TMPDIR+"/ctrl_pts.txt")) syscall("rm "+$TMPDIR+"/ctrl_pts.txt")
if(fexists($TMPDIR+"/project_crism_image.tif1")) syscall("rm "+$TMPDIR+"/project_crism_image.tif1")
if(fexists($TMPDIR+"/project_crism_image.tif2")) syscall("rm "+$TMPDIR+"/project_crism_image.tif2")
if (manual==0) {
for(i=5;i<=dlatlon[1]+1;i+=dlatlon[1]/10) {
for(j=5;j<=dlatlon[2]+1;j+=dlatlon[2]/10){
if(i>dlatlon[1]) {
i=dlatlon[1]
}
if(j>dlatlon[2]) {
j=dlatlon[2]
}
fprintf($TMPDIR+"/ctrl_pts.txt","-gcp %i %i %f %f\n",i,j,latlon[i,j,2],latlon[i,j,1])
}
}
} else {
for (i=1; i<=dlatlon[2]; i+=1) {
if (latlon[1,i]!=0){
fprintf($TMPDIR+"/ctrl_pts.txt","-gcp %i %i %f %f\n",latlon[1,i],latlon[2,i],latlon[3,i],latlon[4,i])
}
}
}
if(HasValue(clat)) {
centerlat=clat
} else {
centerlat=0
}
if(HasValue(clon)) {
centerlon=clon
} else {
if (manual==0) {
centerlon=latlon[dlatlon[1]/2,dlatlon[2]/2,2]
} else {
centerlon=avg(latlon[3])
}
}
if(planet==1) {
#set the projection
system($DV_GDAL_PATH+"/gdal_translate -a_srs GEOGCS[\"Mars 2000\",DATUM[\"D_Mars_2000\",SPHEROID[\"Mars_2000_IAU_IAG\",3396190.0,169.89444722361179]],PRIMEM[\"Greenwich\",0]] `cat "+$TMPDIR+"/ctrl_pts.txt` "+$TMPDIR+"/project_crism_temp_image.tif "+$TMPDIR+"/project_crism_image.tif1")
#project the image using specified projection parameters
#simp
if(proj=="simp") {
syscall($DV_GDAL_PATH+"/gdalwarp -r bilinear "+sampopt+" -srcnodata "+ignore+" -t_srs \"+proj=eqc +lat_ts="+centerlat+" +lat_0="+centerlat+" +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs\" -tps "+$TMPDIR+"/project_crism_image.tif1 "+$TMPDIR+"/project_crism_image.tif2")
}
#sinu
if(proj=="sinu") {
syscall($DV_GDAL_PATH+"/gdalwarp -q -rb "+sampopt+" -srcnodata "+ignore+" -t_srs \"+proj=sinu +lon_0="+centerlon+" +x_0=0 +y_0=0 +a=3396190 +b=3376200 +units=m +no_defs\" -tps "+$TMPDIR+"/project_crism_image.tif1 "+$TMPDIR+"/project_crism_image.tif2")
}
}
if(planet==2) {
syscall($DV_GDAL_PATH+"/gdal_translate -strict -a_srs epsg:4326 `cat "+$TMPDIR+"/ctrl_pts.txt` "+$TMPDIR+"/project_crism_temp_image.tif "+$TMPDIR+"/project_crism_image.tif1")
#simp
if(proj=="simp") {
syscall($DV_GDAL_PATH+"/gdalwarp -r bilinear -tps "+sampopt+" -srcnodata "+ignore+" -t_srs \"+proj=eqc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0="+centerlon+" +units=m +x_0=0 +y_0=0\" -tps "+$TMPDIR+"/project_crism_image.tif1 "+$TMPDIR+"/project_crism_image.tif2")
}
#sinu
if(proj=="sinu") {
syscall($DV_GDAL_PATH+"/gdalwarp -r bilinear "+sampopt+" -srcnodata "+ignore+" -t_srs \"+proj=sinu +lon_0="+centerlon+" +units=m +x_0=0 +y_0=0\" -tps "+$TMPDIR+"/project_crism_image.tif1 "+$TMPDIR+"/project_crism_image.tif2")
}
}
if(planet==3) {
#set the projection
#project the image using specified projection parameters
system($DV_GDAL_PATH+"/gdal_translate -a_srs GEOGCS[\"Moon 2000\",DATUM[\"D_Moon_2000\",SPHEROID[\"Moon_2000_IAU_IAG\",1737400.0,0.0]],PRIMEM[\"Greenwich\",0],UNIT[\"Decimal_Degree\",0.0174532925199433]] `cat "+$TMPDIR+"/ctrl_pts.txt` "+$TMPDIR+"/project_crism_temp_image.tif "+$TMPDIR+"/project_crism_image.tif1")
#project the image using specified projection parameters
#simp
if(proj=="simp") {
syscall($DV_GDAL_PATH+"/gdalwarp -r bilinear "+sampopt+" -srcnodata "+ignore+" -t_srs \"+proj=eqc +lat_ts="+centerlat+" +lat_0="+centerlat+" +x_0=0 +y_0=0 +a=1737400 +b=1737400 +units=m +no_defs\" -tps "+$TMPDIR+"/project_crism_image.tif1 "+$TMPDIR+"/project_crism_image.tif2")
}
#sinu
if(proj=="sinu") {
syscall($DV_GDAL_PATH+"/gdalwarp -q -rb "+sampopt+" -srcnodata "+ignore+" -t_srs \"+proj=sinu +lon_0="+centerlon+" +x_0=0 +y_0=0 +a=1737400 +b=1737400 +units=m +no_defs\" -tps "+$TMPDIR+"/project_crism_image.tif1 "+$TMPDIR+"/project_crism_image.tif2")
}
}
system("convert "+$TMPDIR+"/project_crism_image.tif2 "+$TMPDIR+"/project_crism.pgm")
if(bnum==1) {
outimage=read($TMPDIR+"/project_crism.pgm")
} else {
outimage=cat(outimage,read($TMPDIR+"/project_crism.pgm"),axis=z)
}
}
#syscall("rm "+$TMPDIR+"/ctrl_pts.txt "+$TMPDIR+"/project_crism*")
if (outstruct==1) {
out=get_proj_info()
out.image=outimage
return(out)
} else {
return(outimage)
}
}
define project_crism_index(proj,planet,index,colors,ignore,struct,outstruct,rsstretch) {
#tested cse 3/30/09
if($ARGC!=3 && HasValue(struct)==0) {
printf("\n")
printf("Projects and colorizes CRISM index images\n")
printf("$1 = index image\n")
printf("$2 = data from which the index was created\n")
printf("$3 = ddr for project_crism()\n")
printf("Or\n")
printf("struct = standard CRISM structure with all indices present (from process_crism())\n\n")
printf("Options:\n")
printf("\tproj=projection type (\"sinu\",\"simp\", Default=\"sim\"\n")
printf("\tplanet=use the Mars ellipsoid reference (Default=1)\n")
printf("\tindex=use if more than one index is supplied (Default=1)\n")
printf("\t\te.g. index=5 is the 5th index in the index array in the z direction\n")
printf("\tcolors=color options from colorize() (Default=5)\n")
printf("\tignore=ignore value of the image (Default=0)\n")
printf("\toutstruct=run gdalinfo and return geometry info, Default is 1\n")
printf("\trsstretch=run rsstretch rather than sstretch, Default is 0\n")
printf("NOTE: requires gdalwarp and gdal_translate >1.5.1, specified by $DV_GDAL_PATH\n\n")
printf("modified by c.edwards, created by l.mayorga 3/09\n\n")
return (null)
}
if(HasValue(proj)==0) proj="sinu"
if(HasValue(planet)==0) planet=1
if(HasValue(index)==0) index=1
if(HasValue(colors)==0) colors=5
if(HasValue(ignore)==0) ignore=0
if(HasValue(outstruct)==0) outstruct=1
if(HasValue(rsstretch)==0) rsstretch=0
#parse the inputs
if((HasValue(struct)==1) && type(struct)=="STRUCT") {
indeximg=struct.long.index[,,index]
data=struct.long.data
ddr=struct.long.anc.ddr[,,4:5]
} else if ($ARGC==3) {
indeximg=$1
indeximg=indeximg[,,index]
data=$2
ddr=$3
} else {
printf("No valid data was supplied\n")
return(null)
}
#stretch the specified index
if (rsstretch==1) {
tmp=rsstretch(indeximg,ignore=ignore)
tmp [where (tmp==0 && data[,,50]!=ignore)]=1
} else {
tmp=sstretch(indeximg,ignore=ignore)
tmp [where (tmp==0 && data[,,50]!=ignore)]=1
}
#colorize the index
tmpc=colorize(tmp,ignore=ignore,colors=colors)
#set up the background image
if (dim(data)[3]==438){
bg=avg(data[,,350:420],axis=z)
} else if (dim(data)[3]==55){
bg=avg(data[,,35:45],axis=z)
} else {
return(null)
}
#process and stretch the background image
if (rsstretch==1) {
bg=hpf(bg,size=200,ignore=ignore)
bg=rsstretch(bg,ignore=ignore)
} else {
bg=hpf(bg,size=200,ignore=ignore)
bg=sstretch(bg,ignore=ignore)
}
#convert to hsv and do the overlay
tmpch=rgb2hsv(tmpc)
tmpch[,,3]=(bg+40)/255.
rgb=byte(hsv2rgb(tmpch,maxval=255))
rgb[where (data[,,50]==ignore)]=0
#project the index
prgb=project_crism(ddr,rgb,planet=planet,proj=proj,outstruct=outstruct)
return(prgb)
}
define create_indices(data) {
#tested cse 3/30/09
if(HasValue(data)==0 && $ARGC==1) {
data=$1
} else if (HasValue(data)==0 && $ARGC==0) {
printf("Create a variety of standard indexes for long CRISM FRT/HRL/MSP products\n")
printf("$1=standard CRISM data array (e.g. out.long.data)\n")
printf("Or\n")
printf("data=standard CRISM data array (e.g. out.long.data)\n\n")
printf("Options:\n")
printf("modified by c.edwards, created by j.bandfield, 3/09\n")
printf("now works with short wavelength array data, 10/11\n")
return(null)
}
if(dim(data)[3]==438) {
r1050=avg(data[,,423:438],axis=z)
r1210=avg(data[,,399:414],axis=z)
r1330=avg(data[,,381:396],axis=z)
r1470=avg(data[,,359:374],axis=z)
r1695=avg(data[,,324:340],axis=z)
r1815=avg(data[,,307:322],axis=z)
#r2067=avg(data[,,270:285],axis=z)
#r2530=avg(data[,,198:214],axis=z)
r2100=avg(data[,,270:272],axis=z)
r2400=avg(data[,,224:227],axis=z)
r2210=avg(data[,,253:255],axis=z)
r2120=avg(data[,,267:269],axis=z)
r2140=avg(data[,,264:266],axis=z)
r2250=avg(data[,,248:250],axis=z)
r2250w=avg(data[,,243:257],axis=z)
r2140w=avg(data[,,262:271],axis=z)
r2350w=avg(data[,,230:240],axis=z)
r2290=avg(data[,,241:243],axis=z)
r2350=avg(data[,,232:234],axis=z)
r1930=avg(data[,,296:298],axis=z)
r1985=avg(data[,,287:289],axis=z)
r1857=avg(data[,,307:309],axis=z)
r2067=avg(data[,,275:277],axis=z)
r2530=avg(data[,,205:207],axis=z)
r2330=avg(data[,,235:237],axis=z)
r2230=avg(data[,,250:252],axis=z)
r2390=avg(data[,,226:228],axis=z)
r2600=avg(data[,,194:196],axis=z)
slope=(r1815-r2530)/(2530-1815)
cr2290_cr2330=avg(data[,,236:242],axis=z)/((2310-1815)*slope+r1815)
cr2140_cr2210=avg(data[,,254:265],axis=z)/((2175-1815)*slope+r1815)
r3835=avg(data[,,15:25],axis=z)
r3615=avg(data[,,45:55],axis=z)
# Olivine Index2 from CAT - less slope dependant
r1054=data[,,430]
r1211=data[,,406]
r1329=data[,,388]
r1474=data[,,366]
slope2=clone(slope*0.,z=2)
for(i=1;i<=dim(data)[1];i+=1) {
for(j=1;j<=dim(data)[2];j+=1) {
tmp=fit(data[i,j,225:324],x=create(1,1,100,start=225),type=linear)
slope2[i,j,1]=tmp[1]
slope2[i,j,2]=tmp[2]
}
}
cr1054=430*slope2[,,2]+slope2[,,1]
cr1211=406*slope2[,,2]+slope2[,,1]
cr1329=388*slope2[,,2]+slope2[,,1]
cr1474=366*slope2[,,2]+slope2[,,1]
indexol2=(((cr1054-r1054)/cr1054)*0.1) + (((cr1211-r1211)/cr1211)*0.1) + (((cr1329-r1329)/cr1329)*0.4)+(((cr1474-r1474)/cr1474)*0.4)
# Olivine Index
index=r1695/(r1050*0.1+r1210*0.1+r1330*0.4+r1470*0.4)-1
# LCP Index
index=cat(index,(r1330-r1050)/(r1330+r1050)*(r1330-r1815)/(r1330+r1815),axis=z)
# HCP Index
index=cat(index,(r1470-r1050)/(r1470+r1050)*(r1470-r2067)/(r1470+r2067),axis=z)
# Ferric Coating
index=cat(index,(r1815-r2530)/(2530-1815),axis=z)
# BD2290 MG,Fe-OH minerals
index=cat(index,1-(r2290/(0.6*r2250+0.4*r2350)),axis=z)
# D2300 hydrated minerals - particularly phyllosilicates
index=cat(index,(1-(cr2290_cr2330-cr2140_cr2210)),axis=z)
# BD2210 Al-OH minerals
index=cat(index,1-(r2210/(((4./11.)*r2140)+((7./11.)*r2250))),axis=z)
# BD1900 H20
# BD1900=1.-((.5*(r1930+r1985))/(((167./210.)*r1857)+((43./210.)*r2067)))
# Alternate from CAT code
index=cat(index,1.- (r1985 + r1930) / (r2067 + r1857),axis=z)
# SINDEX Water containing minerals
index=cat(index,1-((r2100+r2400)/(2*r2290)),axis=z)
# BD2100 Monohydrated sulfates
index=cat(index,1-((r2120+r2140)*0.5)/(0.375*r1930+0.625*r2250),axis=z)
# BDCARB Carbonate Index
index=cat(index,1-((r2330/(0.375*r2230+0.625*r2390))*(r2530/(0.33*r2390+0.67*r2600)))^0.5,z)
# BD3000 H20 3
index=cat(index,1-(r2530/(0.33*r2390+0.67*r2600)),z)
index=cat(index,1-(r2330/(0.375*r2230+0.625*r2390)),z)
index=cat(index,1-(r3615/(r3835)),z)
# Hydrated Silica Index
index=cat(index,r2140w/r2250w,z)
# Zeolite Index
zeo=avg(data[,,245:255],axis=z)-avg(data[,,200:210],axis=z)
index=cat(index,zeo,z)
index=cat(index,indexol2,z)
} else if(dim(data)[3]==55) {
r1050=avg(data[,,53:55],axis=z)
r1210=avg(data[,,49:52],axis=z)
r1330=avg(data[,,45:48],axis=z)
r1470=avg(data[,,40:43],axis=z)
r1695=avg(data[,,35],axis=z)
r1815=avg(data[,,34],axis=z)
#r2067=avg(data[,,28:29],axis=z)
#r2530=avg(data[,,14],axis=z)
r2400=avg(data[,,17],axis=z)
r2210=avg(data[,,24],axis=z)
r2120=avg(data[,,27],axis=z)
r2140=avg(data[,,26],axis=z)
r2250=avg(data[,,22:23],axis=z)
r2250w=avg(data[,,22:24],axis=z)
r2140w=avg(data[,,26:27],axis=z)
r2350w=avg(data[,,18:20],axis=z)
r2290=avg(data[,,21],axis=z)
r2350=avg(data[,,18],axis=z)
r1930=avg(data[,,32],axis=z)
r1985=avg(data[,,30:31],axis=z)
r2067=avg(data[,,28],axis=z)
r2530=avg(data[,,14],axis=z)
r2330=avg(data[,,19],axis=z)
r2230=avg(data[,,23],axis=z)
r2390=avg(data[,,17],axis=z)
r2600=avg(data[,,13],axis=z)
slope=(r1815-r2530)/(2530-1815)
cr2290_cr2330=avg(data[,,19:21],axis=z)/((2310-1815)*slope+r1815)
cr2140_cr2210=avg(data[,,24:25],axis=z)/((2175-1815)*slope+r1815)
r3615=avg(data[,,4],axis=z)
# Olivine Index
index=r1695/(r1050*0.1+r1210*0.1+r1330*0.4+r1470*0.4)-1
index [where (r1695==0)]=0
# LCP Index
index=cat(index,(r1330-r1050)/(r1330+r1050)*(r1330-r1815)/(r1330+r1815),axis=z)
# HCP Index
index=cat(index,(r1470-r1050)/(r1470+r1050)*(r1470-r2067)/(r1470+r2067),axis=z)
# Ferric Coating
index=cat(index,(r1815-r2530)/(2530-1815),axis=z)
# BD2290 MG,Fe-OH minerals
index=cat(index,1-(r2290/(0.6*r2250+0.4*r2350)),axis=z)
# D2300 hydrated minerals - particularly phyllosilicates
index=cat(index,(1-(cr2290_cr2330-cr2140_cr2210)),axis=z)
# BD2210 Al-OH minerals
index=cat(index,1-(r2210/(((4./11.)*r2140)+((7./11.)*r2250))),axis=z)
# BD1900 H20
# BD1900=1.-((.5*(r1930+r1985))/(((167./210.)*r1857)+((43./210.)*r2067)))
# Alternate from CAT code
# index=cat(index,1.- (r1985 + r1930) / (r2067 + r1857),axis=z)
# SINDEX Water containing minerals
# index=cat(index,1-((r2100+r2400)/(2*r2290)),axis=z)
# BD2100 Monohydrated sulfates
index=cat(index,1-((r2120+r2140)*0.5)/(0.375*r1930+0.625*r2250),axis=z)
# BDCARB Carbonate Index
index=cat(index,1-((r2330/(0.375*r2230+0.625*r2390))*(r2530/(0.33*r2390+0.67*r2600)))^0.5,z)
index=cat(index,1-(r2530/(0.33*r2390+0.67*r2600)),z)
index=cat(index,1-(r2330/(0.375*r2230+0.625*r2390)),z)
# index=cat(index,1-(r3615/(r3835)),z)
# Hydrated Silica Index
index=cat(index,r2140w/r2250w,z)
# Zeolite Index
zeo=avg(data[,,22:24],axis=z)-avg(data[,,14],axis=z)
index=cat(index,zeo,z)
# If the short wavelength array is present, then create the short indices
} else if(dim(data)[3]==107) {
index=data[,,36]
index=cat(index,data[,,27],z)
index=cat(index,data[,,21],z)
r440=data[,,13]
r530=data[,,27]
r600=data[,,38]
r709=data[,,54]
r800=data[,,68]
r920=data[,,87]
r984=data[,,97]
#BD530 Ferric Minerals
#1 − (R530/(a*R709 + b*R440))
#BDC = 1 − RC/(a*RS + b*RL), where a = 1 − b and b = (λC − λS)/(λL − λS).
b=(530-440)/(709-440)
a=1-b
index=cat(index,(1-r530/(a*r709+b*r440)),z)
#SH600 Coatings
#R600/(a*R530 + b*R709)
b=(600-530)/(709-530)
a=1-b
index=cat(index,r600/(a*r530 + b*r709),z)
#BD1000 Iron Minerals (using BD920 as a simple to calculate substitute)
#1 - ( cube(*,*,R920) / (b * cube(*,*,R800) + a * cube(*,*,R984) ) )
b=(920-800)/(984-800)
a=1-b
index=cat(index,1-(r920/(a*r800 + b*r984)),z)
# Short wavelength multispectral data (closest wavelengths available)
} else if(dim(data)[3]==19) {
index=data[,,5]
index=cat(index,data[,,4],z)
index=cat(index,data[,,3],z)
r440=data[,,3]
r530=data[,,4]
r600=data[,,5]
r709=data[,,8]
r800=data[,,11]
r920=data[,,15]
r984=data[,,17]
#BD530 Ferric Minerals
#1 − (R530/(a*R709 + b*R440))
#BDC = 1 − RC/(a*RS + b*RL), where a = 1 − b and b = (λC − λS)/(λL − λS).
b=(530-440)/(709-440)
a=1-b
index=cat(index,(1-r530/(a*r709+b*r440)),z)
#SH600 Coatings
#R600/(a*R530 + b*R709)
b=(600-530)/(709-530)
a=1-b
index=cat(index,r600/(a*r530 + b*r709),z)
#BD1000 Iron Minerals (using BD920 as a simple to calculate substitute)
#1 - ( cube(*,*,R920) / (b * cube(*,*,R800) + a * cube(*,*,R984) ) )
b=(920-800)/(984-800)
a=1-b
index=cat(index,1-(r920/(a*r800 + b*r984)),z)
}
return(float(index))
}
define create_indices_key() {
#tested cse 3/30/09
if($ARGC!=1) {
printf("Create the index key based for the FRT/HRL/MSP indexes\n")
printf("$1=array of indexes\n")
printf("modified by c.edwards, created by j.bandfield, 3/09\n")
return(null)
}
indices=$1
if (dim(indices)[3]>=16){
index_key=cat("Olivine Index","LCP Index","HCP Index","Ferric Coating","BD2290 MG, Fe-OH Minerals","D2300 Hydrates Minerals - Particularily Phyllosilicates","BD2210 AL-OH Minerals","BD1900 H20 - Alternate from CAT code","SINDEX Water Containing Minerals","BD2100 Monohydrated Sulfates","BDCARB Carbonate Index","BD3000 H20 3","Unknown1","Unknown2","Hydrated Silica Index","Zeolite Index","Olivine Index2",y)
} else if (dim(indices)[3]==13){
index_key=cat("Olivine Index","LCP Index","HCP Index","Ferric Coating","BD2290 MG, Fe-OH Minerals","D2300 Hydrates Minerals - Particularily Phyllosilicates","BD2210 AL-OH Minerals","BD2100 Monohydrated Sulfates","BDCARB Carbonate Index","BD3000 H20 3","Unknown","Hydrated Silica Index","Zeolite Index",y)
} else if (dim(indices)[3]==6){
index_key=cat("Red","Green","Blue","BD530 Ferric Minerals","SH600 Coatings","BD920 Fe-Mineralogy",y)
} else {
return(null)
}
return(index_key)
}
define crism_speclib(spec_lib,srch) {
#tested-cse 3/30/09
if($ARGC!=1){
printf("\nSearch through the standard VNIR spectral library\n")
printf("$1 = whatever you are searching for, usually in all caps, STRING\n")
printf("srch = location to search in (default=lib.class)\n")
printf("spec_lib = variable the spectral library is stored (default will be read from $DV_SCRIPT_FILES/crism_speclib.hdf\n")
printf("Ex:\n")
printf("a=crism_speclib(\"OLIVINE\")\n")
printf("a=crism_speclib(\"OLIVINE\",srch=\"lib.name\")\n")
printf("modified by c.edwards, created by l.mayorga 3/09")
printf("\n\n")
return(null)
}
in=$1
if(HasValue(spec_lib)==0) {
verbose=0
spec_lib=$DV_SCRIPT_FILES+"/crism_speclib.hdf"
lib=read(spec_lib)
verbose=1
}
if(HasValue(srch)==0) srch="lib.class"
printf("searching in %s\n\n",srch)
srch=eval(srch)
out=struct(index,name,data)
out.index=int(0)
out.name="\"\""
out.data=struct()
out.id="\"\""
out.type="\"\""
out.version_id="\"\""
out.specimen_name="\"\""
out.specimen_description="\"\""
out.min_particle_size=float(0)
out.max_particle_size=float(0)
out.collection_location="\"\""
out.current_location="\"\""
out.reference="\"\""
counter=1
inlen=strlen(in)
nlines=dim(strlen(srch))[2]
nchars=strlen(srch)
for (i=1; i<=nlines; i+=1) {
for (j=1; j<=int((nchars[,i]-inlen)); j+=1) {
if (srch[j:(j+(inlen-1)),i]==in) {
out.index=out.index//int(i)
out.name=cat(out.name,lib.class[,i],y)
out.id=cat(out.id,lib.id[,i],y)
out.type=cat(out.type,lib.type[,i],y)
out.version_id=cat(out.version_id,lib.version_id[,i],y)
out.specimen_name=cat(out.specimen_name,lib.specimen_name[,i],y)
out.specimen_description=cat(out.specimen_description,lib.specimen_description[,i],y)
out.min_particle_size=cat(out.min_particle_size,lib.min_particle_size[,i],y)
out.max_particle_size=cat(out.max_particle_size,lib.max_particle_size[,i],y)
out.collection_location=cat(out.collection_location,lib.collection_location[,i],y)
out.current_location=cat(out.current_location,lib.current_location[,i],y)
out.reference=cat(out.reference,lib.reference[,i],y)
add_struct(object=out.data, name=sprintf("id%.4d_%d",i,counter),value=lib.data[i])
j=nchars[,i]
counter=counter+1
}
}
}
if (dim(out.index)[1]==1) {
return(0)
} else {
out.index=out.index[2:]
out.name=out.name[,2:]
out.id=out.id[,2:]
out.type=out.type[,2:]
out.version_id=out.version_id[,2:]
out.specimen_name=out.specimen_name[,2:]
out.specimen_description=out.specimen_description[,2:]
out.min_particle_size=out.min_particle_size[,2:]
out.max_particle_size=out.max_particle_size[,2:]
out.collection_location=out.collection_location[,2:]
out.current_location=out.current_location[,2:]
out.reference=out.reference[,2:]
return(out)
}
}
define get_crism_list(maxlat,minlat,maxlon,minlon,targeted) {
#tested-cse 3/30/09
if(HasValue(maxlat)==0 || HasValue(minlat)==0 || HasValue(maxlon)==0 || HasValue(minlon)==0) {
printf("Returns full resolution targets for a defined lat/lon region\n")
printf("Required:\n")
printf("\tmaxlat=maximum latitude\n")
printf("\tminlat=minimum latitude\n")
printf("\tmaxlon=maximum longitude\n")
printf("\tminlon=minimum longitude\n")
printf("Optional:\n")
printf("\ttargeted=selected targeted images (Default=1)\n\n")
printf("j.bandfield 10/08\n")
return(null)
}
mysql_login="mysql -s -D crism -h octuco.ess.washington.edu -P "+$DV_CRISM_DB_PORT+" -u "+$DV_CRISM_DB_USER+" -p"+$DV_CRISM_DB_PASS+" -e"
if (maxlon >180) {
maxlon=maxlon-360
minlon=minlon-360
}
if (HasValue(targeted)==0) targeted=1
# Only full res targets for now
# Get filename list
if (targeted==0) {
list=syscall(mysql_login+" 'select ddr.id,averaging_width,columns,rows,ullat,urlat,lllat,lrlat,ullon,urlon,lllon,lrlon from trdr,ddr where trdr.id=ddr.id and trdr.sensor_id=ddr.sensor_id and trdr.path like \"\%IF\%\" and (trdr.path like \"\%MSS\%\" or trdr.path like \"\%MSP\%\" or trdr.path like \"\%MSW\%\") and lllat<"+maxlat+" and ullat>"+minlat+" and urlon<"+maxlon+" and ullon>"+minlon+" and trdr.sensor_id like \"L\";'")
}
if (targeted!=0) {
list=syscall(mysql_login+" 'select ddr.id,averaging_width,columns,rows,ullat,urlat,lllat,lrlat,ullon,urlon,lllon,lrlon from trdr,ddr where trdr.id=ddr.id and trdr.sensor_id=ddr.sensor_id and trdr.path like \"\%IF\%\" and (trdr.path like \"\%FRT\%\" or trdr.path like \"\%HRL\%\" or trdr.path like \"\%HRS\%\") and lllat<"+maxlat+" and ullat>"+minlat+" and urlon<"+maxlon+" and ullon>"+minlon+" and trdr.sensor_id like \"L\" and averaging_width<=2;'")
}
return(list)
}
define get_proj_info() {
out=struct()
out.text=syscall($DV_GDAL_PATH+"/gdalinfo "+$TMPDIR+"/project_crism_image.tif2")
nlines=dim(strlen(out.text))[2]
i=1
while (i<=nlines) {
if (out.text[5:14,i]=="PROJECTION") {
out.projection=out.text[5:,i]
}
if (out.text[1:10,i]=="Pixel Size") {
out.sampling=atof(out.text[15:25,i])
}
if (out.text[16:34,i]=="longitude of center") {
out.clon=atof(out.text[37:45,i])
}
if (out.text[1:6,i]=="Corner") {
cline=i+1
}
if (out.text[1:6,i]=="Origin") {
originline=out.text[11:,i]
}
i+=1
}
chars=strlen(originline)-1
originline=originline[1:chars]
i=1
commatrig=0
while (i<=chars) {
if (originline[i]==",") {
out.offset=atod(originline[1:(i-1)])//atod(originline[(i+1):])
out.sampleoffset=int(atod(originline[1:(i-1)])/out.sampling)
out.lineoffset=-int(atod(originline[(i+1):])/out.sampling)
}
i+=1
}
out.lon=create(1,5,1,start=0,step=0,format=float)
out.lat=out.lon
for (j=0; j<=4; j+=1) {
chars=strlen(out.text[,cline+j])
i=1
commatrig=0
d1start=0
m1start=0
s1start=0
s1stop=0
d2start=0
m2start=0
s2start=0
while (i<=chars) {
if (out.text[i:i+2,cline+j]==") (") {
d1start=i+3
}
if (out.text[i,cline+j]=="d" && d1start!=0 && commatrig==0) {
d1stop=i-1
m1start=i+1
}
if (out.text[i,cline+j]=="'" && m1start!=0 && commatrig==0) {
m1stop=i-1
s1start=i+1
}
if (out.text[i,cline+j]=="\"" && s1start!=0 && commatrig==0) {
s1stop=i-1
if (out.text[i+1,cline+j]=="E") {
east=1
} else {
east=-1
}
}
if (out.text[i,cline+j]=="," && s1stop!=0) {
d2start=i+1
commatrig=1
}
if (out.text[i,cline+j]=="d" && d2start!=0 && commatrig==1) {
d2stop=i-1
m2start=i+1
}
if (out.text[i,cline+j]=="'" && m2start!=0 && commatrig==1) {
m2stop=i-1
s2start=i+1
}
if (out.text[i,cline+j]=="\"" && s2start!=0 && commatrig==1) {
s2stop=i-1
if (out.text[i+1,cline+j]=="N") {
north=1
} else {
north=-1
}
}
i+=1
}
deg1=atof(out.text[d1start:d1stop,cline+j])
deg2=atof(out.text[d2start:d2stop,cline+j])
min1=atof(out.text[m1start:m1stop,cline+j])
min2=atof(out.text[m2start:m2stop,cline+j])
sec1=atof(out.text[s1start:s1stop,cline+j])
sec2=atof(out.text[s2start:s2stop,cline+j])
out.lon[,j+1]=deg1+min1/60.+sec1/3600.
out.lat[,j+1]=deg2+min2/60.+sec2/3600.
out.lat[,j+1]=out.lat[,j+1]*north
if (east==-1) {
out.lon[,j+1]=360-out.lon[,j+1]
}
}
minlat=min(out.lat)
minlon=min(out.lon)
maxlat=max(out.lat)
maxlon=max(out.lon)
out.ul=maxlat//minlon
out.ur=maxlat//maxlon
out.ll=minlat//minlon
out.lr=minlat//maxlon
return(out)
}
define get_crism_old(ra,raw,sl,force,verb,ignore,set_ignore) {
#tested-cse 3/30/09
if($ARGC!=1) {
printf("Download any CRISM image, short and long, radiance or I/F\n\n")
printf("$1=product id (truncated from full ID)\n\n")
printf("Optional:\n")
printf("\tsl=process both the short and long data product (Default=0, long only)\n")
printf("\tra=download the radiance product instead of I/F (Default=0)\n")
printf("\traw=return just the output structure (other options will not work with this, Default=0)\n")
printf("\tforce=force the redownload of products (Default=0), uses $DV_CRISM_CACHE\n")
printf("\tverb=set the verbosity (Default=0)\n")
printf("\tignore=the ignore value (Default=0)\n\n")
printf("\tset_ignore=set the value of value of the image (Default=1)\n\n")
printf("NOTE: Also uses $DV_CRISM_DB_PORT, $DV_CRISM_DB_USER, and optionally $DV_CRISM_DB_PASS\n\n")
printf("modified by c.edwards, adapted from j.bandfield 3/09\n")
return(null)
}
product=$1
global(verbose)
verbold=verbose
printf("Running get CRISM...")
if (HasValue(ra)==0) ra=0
if (HasValue(raw)==0) raw=0
if (HasValue(sl)==0) sl=0
if (HasValue(verb)==0) verbose=0
if (HasValue(force)==0) force=1
if (HasValue(ignore)==0) ignore=0
if (HasValue(set_ignore)==0) set_ignore=0
mysql_login="mysql -s -D crism -h octuco.ess.washington.edu -P "+$DV_CRISM_DB_PORT+" -u "+$DV_CRISM_DB_USER+" -p"+$DV_CRISM_DB_PASS+" -e"
if (ra==0) {
files=syscall(mysql_login+" 'select ddr.volume_id,ddr.path,trdr.volume_id, trdr.path from trdr,ddr where ddr.id=\""+product+"\" and trdr.id=ddr.id and trdr.sensor_id=ddr.sensor_id and trdr.path like \"\%IF\%\" group by trdr.path;'")
} else {
files=syscall(login+" 'select ddr.volume_id,ddr.path,trdr.volume_id,trdr.path from trdr,ddr where ddr.id=\""+product+"\" and trdr.id=ddr.id and trdr.sensor_id=ddr.sensor_id and trdr.path like \"\%RA\%\" group by trdr.path;'")
}
if(HasValue(files)==0) {
printf("\nWrong password or incorrect image ID\n")
printf("Please try again\n\n")
verbose=verb_old
return(null)
}
ddrvolnum=files[9:10,1]
volnum=files[81:82,1]
tmpdir=sprintf("%s/%s",$DV_CRISM_CACHE,product)
printf("%s\n",tmpdir)
system(sprintf("mkdir -p %s",tmpdir))
system(sprintf("chmod 777 %s",tmpdir))
#system("cd "+tmpdir+"; wget -nd -nc http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:71,1]+" "+tmpdir)
#system("cd "+tmpdir+"; wget -nd -nc http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:68,1]//"IMG "+tmpdir)
#get_file("http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:71,1],tmpdir,force=force)
#get_file("http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:68,1]//"IMG",tmpdir,force=force)
#system("cd "+tmpdir+"; wget -nd -nc http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:144,1]+" "+tmpdir)
#system("cd "+tmpdir+"; wget -nd -nc http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:141,1]//"IMG "+tmpdir)
#get_file("http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:144,1],tmpdir,force=force)
#get_file("http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:141,1]//"IMG",tmpdir,force=force)
ddrpath1="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-6-ddr-v1/mrocr_10"//ddrvolnum//"/ddr"//"/"//files[16:71,1]
ddrpath3="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-6-ddr-v1/mrocr_10"//ddrvolnum//"/ddr"//"/"//files[16:68,1]//"IMG"
system(sprintf("cd %s; %s",tmpdir,ddrpath1))
system(sprintf("cd %s; %s",tmpdir,ddrpath3))
trrpath1="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism-trdr-1/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:144,1]
trrpath3="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism-trdr-1/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:141,1]//"IMG"
system(sprintf("cd %s; %s",tmpdir,trrpath1))
system(sprintf("cd %s; %s",tmpdir,trrpath3))
if (sl!=0) {
system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:71,2]+" "+tmpdir)
system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:68,2]//"IMG "+tmpdir)
#get_file("http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:71,2],tmpdir,force=force)
#get_file("http://pds-geosciences.wustl.edu/crism-ddr01/mro-m-crism-6-ddr-v1/mrocr_1001/ddr"//"/"//files[16:68,2]//"IMG",tmpdir,force=force)
system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:144,2]+" "+tmpdir)
system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:141,2]//"IMG "+tmpdir)
#get_file("http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:144,2],tmpdir,force=force)
#get_file("http://pds-geosciences.wustl.edu/crism-trdr"//volnum//"/mro-m-crism-3-rdr-targeted-v1/mrocr_20"//volnum//"/trdr"//"/"//files[89:141,2]//"IMG",tmpdir,force=force)
}
if (sl!=0) {
data=struct(short,long,ddrs,ddrl)
} else {
data=struct(long,ddrl)
}
list=syscall(sprintf("ls %s/*%s*.LBL",tmpdir,product))
if (sl!=0) {
shortl=syscall("nl "+list[,4]+" | grep HKP | awk '{print $1-3}'")
if(HasValue(shortl)==1) {
system("head -n "+shortl[,1]+" "+list[,4]+" > "+tmpdir+"/temp")
system("cat "+tmpdir+"/temp "+tmpdir+"/temp2> "+list[,4])
}
data.ddrs=load_pds(list[,2])
data.short=load_pds(list[,4])
ll=syscall("nl "+list[,3]+" | grep HKP | awk '{print $1-3}'")
if(HasValue(ll)==1) {
system("head -n "+ll[,1]+" "+list[,3]+" > "+tmpdir+"/temp")
system("echo END > "+tmpdir+"/temp2")
system("cat "+tmpdir+"/temp "+tmpdir+"/temp2> "+list[,3])
}
data.ddrl=load_pds(list[,1])
data.long=load_pds(list[,3])
}
if (sl==0) {
list=grep(list,"L_")
ll=syscall("nl "+list[,2]+" | grep HKP | grep L_ | awk '{print $1-3}'")
if(HasValue(ll)==1) {
system("head -n "+ll[,1]+" "+list[,2]+" > "+tmpdir+"/temp")
system("echo END > "+tmpdir+"/temp2")
system("cat "+tmpdir+"/temp "+tmpdir+"/temp2> "+list[,2])
}
data.ddrl=load_pds(list[,1])
data.long=load_pds(list[,2])
}
#the way it should be done
# if (sl!=0) {
# data.ddrl=load_pds(list[,1])
# data.long=load_pds(list[,3])
# } else {
# data.ddrl=load_pds(list[,1])
# data.long=load_pds(list[,2],data=0)
# data.long.file.image=read(list[,2])
# }
# Get wavelength arrays
wl_filename1=get_struct(data.long,"mro\:wavelength_file_name")
wl_len1=strlen(wl_filename1)
wl_filename1=wl_filename1[2:(wl_len1-4)]
#system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename1//"LBL "+tmpdir)
#system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename1//"IMG "+tmpdir)
# get_file("http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename1//"LBL",tmpdir,force=force)
# get_file("http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename1//"IMG",tmpdir,force=force)
wl1="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename1//"LBL"
wl3="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename1//"IMG"
system(sprintf("cd %s; %s",tmpdir,wl1))
system(sprintf("cd %s; %s",tmpdir,wl3))
list=syscall(sprintf("ls %s/CDR*WA*L_*.LBL",tmpdir))
data.long.wl_array=load_pds(list[,1])
if (sl!=0) {
wl_filename2=get_struct(data.short,"mro\:wavelength_file_name")
wl_len2=strlen(wl_filename2)
wl_filename2=wl_filename2[2:(wl_len2-4)]
system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename2//"LBL",tmpdir)
system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename2//"IMG",tmpdir)
#get_file("http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename2//"LBL",tmpdir,force=force)
#get_file("http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/wa/"//wl_filename2//"IMG",tmpdir,force=force)
list=syscall(sprintf("ls %s/CDR*WA*S_*.LBL",tmpdir))
data.short.wl_array=load_pds(list[,1])
}
# Get atmospheric transmission arrays
# Use these for the atmospheric correction (long wave only)
at_filename1=get_struct(data.long,"mro\:wavelength_file_name")
at_len1=strlen(at_filename1)
at_filename1=at_filename1[2:(at_len1-4)]
at_filename1=at_filename1[1:5]//"0000000000_AT"//wl_filename1[19:]
#system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//at_filename1//"LBL "+tmpdir)
#system("cd "+tmpdir+"; wget -nc -nd http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//at_filename1//"IMG "+tmpdir)
#get_file("http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//at_filename1//"LBL",tmpdir,force=force)
#get_file("http://pds-geosciences.wustl.edu/crism-edr01/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//at_filename1//"IMG",tmpdir,force=force)
at1="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//at_filename1//"LBL"
at3="wget -nc -nd http://pds-geosciences.wustl.edu/mro-crism/mro-m-crism-2-edr-v1/mrocr_0001/cdr/at/"//at_filename1//"IMG"
system(sprintf("cd %s; %s",tmpdir,at1))
system(sprintf("cd %s; %s",tmpdir,at3))
list=syscall(sprintf("ls %s/CDR*AT*L_*.LBL",tmpdir))
data.long.at_array=load_pds(list[,1])
verbose=verbold
if (raw==1) {
return(data)
}
if (raw!=1) {
out=struct(long)
out.long=struct(anc,wl_array,at_array,id,data)
out.long.anc=struct(ls,solar_dist,ddr,ddrkey)
out.long.anc.summing=data.long.pixel_averaging_width
out.long.data=remove_struct(data.long.file.image,"data")
if(set_ignore==1) {
out.long.data[where out.long.data>100 || out.long.data<-100]=ignore
}
out.long.id=data.long.file.ptr_to_rownum_table[,1]
out.long.wl_array=remove_struct(data.long.wl_array.file.image,"data")
out.long.at_array=remove_struct(data.long.at_array.file.image,"data")
out.long.anc.ls=data.ddrl.solar_longitude
out.long.anc.solar_dist=data.ddrl.solar_distance
out.long.anc.ddrkey=data.ddrl.file.image.band_name
out.long.anc.ddr=remove_struct(data.ddrl.file.image,"data")
if (sl!=0) {
out.short=struct(anc,wl_array,id,data)
out.short.anc=struct(ls,solar_dist,ddr,ddrkey)
out.short.data=remove_struct(data.short.file.image,"data")
if(set_ignore==1) {
out.short.data[where out.short.data>100 || out.short.data<-100]=ignore
}
out.short.id=data.short.file.ptr_to_rownum_table[,1]
out.short.wl_array=remove_struct(data.short.wl_array.file.image,"data")
out.short.anc.ls=data.ddrs.solar_longitude
out.short.anc.solar_dist=data.ddrs.solar_distance
out.short.anc.ddrkey=data.ddrs.file.image.band_name
out.short.anc.ddr=remove_struct(data.ddrs.file.image,"data")
}
}
printf("Done.\n")
return(out)
}