#Original Function Location: /themis/lib/dav_lib/library/themis_atm.dvrc
Source Code for Function: "tes_fp()" in "themis_atm.dvrc" (Beta)
themis_atm_version=1.01
# 1.01 Remove spectra that don't have all 10 bands.
# 1.01 Add option for gap filling for fully automated atm correction method.
# 1.01 Save atmosphere component.
# THEMIS semi-automatic and fully automated atmospheric correction tools THEMIS_atm.dvrc
# themis_tes_overlay -- Search TES footprints over a THEMIS image using specific criteria
# themis_atm_corr -- Fully automated atmospheric attenuation removal for TES masked THEMIS image
# tes_fp -- Get coordinates and topography information of TES data overlapping THEMIS image
# refine_tes_fp -- Find TES pixels used as training data
# themis_rematm -- Traditional THEMIS atmospheric attenuation removal
define themis_tes_overlay(ock_start,ock_end,minlat,maxlat,em_angle,in_angle,temperature,albedo,dust,waterice,auto){
if($ARGC==0){
printf("\nAccepts standard mysql queries on the tes database at the MSFF for a single THEMIS image overlap\n")
printf("Returns the results of the query with constraints below\n")
printf("$1 = THEMIS image ID, for example 'I01001001'\n")
printf("$2 = select items separated by comma from TES Layer options\n")
printf("auto = 0/1, using TES select criteria for semi-automatic (0) or full atmoated correction (1), default=1\n")
printf("Default arguments: contains st_astext, emissivity, ock, ick, detector\n")
printf("\nOptional arguments for TES selection criteria: \n")
printf("ock_start = minimum orbit counter keeper (optional, default is 1583)\n")
printf("ock_end = maximum orbit counter keeper (optional, default is 10000)\n")
printf("minlat = minimum latitude (optional, default is -60)\n")
printf("maxlat = maximum latitude (optional, default is 60)\n")
printf("em_angle = maximum emission angle (optional, default is 10)\n")
printf("in_angle = maximum incidence angle (optional, default is 80)\n")
printf("temperature = minimum target temperature (optional, default is 250)\n")
printf("albedo = maximum albedo (optional, default is 0.25)\n")
printf("dust = maximum atmospheric dust opacities (optional, default is 0.25)\n")
printf("waterice = maximum water ice opacities (optional, default is 0.15)\n")
printf("If the preprocessed cub file in the directory, it will display the TES coverage over THEMIS image\n\n")
return(null)
}
#jmarsFields = "obs_sclk, orbit_number, orbit_counter_keeper, instrument_time_count, mirror_pointing_angle, imc_count, observation_type,"
#+ "scan_length, data_packet_type, schedule_type, spectrometer_gain, visual_bolometer_gain, preprocessor_detector_number,"
#+ "detector_mask, observation_classification, mission_phase, intended_target, tes_sequence, neon_lamp_status, timing_accuracy,"
#+ "classification_value, quality, hga_motion, solar_panel_motion, algor_patch, imc_patch, momentum_desaturation, "
#+ "equalization_table, primary_diagnostic_temperatures, fft_start_index, is_aerobrake, surface_pressure, nadir_temperature_profile,"
#+ "co2_continuum_temp, spectral_surface_temperature, temperature_profile_residual, nadir_opacity, nadir_opacity_residual,"
#+ "co2_downwelling_flux, total_downwelling_flux, temperature_profile_rating, atmospheric_opacity_rating, atmospheric_calibration_id,"
#+ "temporal_integration_scan_number, raw_visual_bolometer, detector, raw_thermal_bolometer, calibrated_visual_bolometer,"
#+ "lambert_albedo, bolometric_thermal_inertia, bolometric_brightness_temp, visual_bol_calibration_id, thermal_bol_calibration_id,"
#+ "bolometric_inertia_rating, bolometer_lamp_anomaly, longitude, longitude_00, latitude, phase, emission, incidence, planetary_phase,"
#+ "heliocentric_lon, sub_sc_lon, sub_sc_lon_00, sub_sc_lat, sub_solar_lon, sub_solar_lon_00, sub_solar_lat, target_distance,"
#+ "height, altitude, local_time, solar_distance, angular_semidiameter, version_id, spectral_mask, compression_mode, detector_temperature,"
#+ "target_temperature, spectral_thermal_inertia, radiance_calibration_id, major_phase_inversion, algor_risk, calibration_failure,"
#+ "calibration_quality, spectrometer_noise, spectral_inertia_rating, detector_mask_problem, raw_radiance, calibrated_radiance,"
#+ "emissivity, perc_bas1, perc_bas2, perc_hem, perc_dust1, perc_dust2, perc_ice1, perc_ice2, perc_bb, rms_error, tot_dust, tot_ice, hem_index,"
#+ "lw_depth, sw_depth, lw_ratio, atm_adj_emiss, mars_year "
verbose=0
id=$1
if($ARGC==1) {
query="st_astext(fp),emissivity,orbit_counter_keeper,instrument_time_count,detector"
} else {
query=$2
query="st_astext(fp),emissivity,orbit_counter_keeper,instrument_time_count,detector,"+query
}
if(HasValue(auto)==0) auto=1
# Default TES criteria for semi-automatic correction
if(auto==0) {
if(HasValue(ock_start)==0) ock_start=1583
if(HasValue(ock_end)==0) ock_end=7000
if(HasValue(minlat)==0) minlat=-60
if(HasValue(maxlat)==0) maxlat=60
if(HasValue(em_angle)==0) em_angle=10
if(HasValue(in_angle)==0) in_angle=80
if(HasValue(temperature)==0) temperature=255
if(HasValue(albedo)==0) albedo=0.2
if(HasValue(dust)==0) dust=0.15
if(HasValue(waterice)==0) waterice=0.08
}
# Default TES criteria for fully automatic correction
if(auto==1) {
if(HasValue(ock_start)==0) ock_start=1583
if(HasValue(ock_end)==0) ock_end=10000
if(HasValue(minlat)==0) minlat=-60
if(HasValue(maxlat)==0) maxlat=60
if(HasValue(em_angle)==0) em_angle=10
if(HasValue(in_angle)==0) in_angle=80
if(HasValue(temperature)==0) temperature=250
if(HasValue(albedo)==0) albedo=0.25
if(HasValue(dust)==0) dust=0.25
if(HasValue(waterice)==0) waterice=0.15
}
# Checks to see if its a valid format in the form of Ixxxxxxxx #
if(length(id)!=9){
printf("invalid image_ID, must be 9 characters long\n\n")
return(null)
}
# Checks to see if the beginning filename starts appropriately #
if(id[1,1]!="I"){
printf("invalid image_ID, must begin with either an I \n\n")
return(null)
}
# Checks to see if image_id given exists in database #
temp_result = read_lines(themis3db("select file_id from thm3_header where file_id='"+id+"';",psql=1,v=0))
if(length(temp_result)==0){
printf("\nImage ID supplied does not exist in database\n\n")
return(null)
}
# Substitute % for %25
queryout=""
while(strstr(query,'%')!=0){
pos=strstr(query,"%")
queryout=queryout//query[:pos-1]+"%25"
query=query[pos+1:]
}
queryout=queryout//query
# Make a call to pgis_outline from thmteam.pgisgeom, no need to use ST_astext here #
url_outline=(themis3db("select pgis_outline from thmteam.pgisgeom where file_id='"+id+"';",psql=1,header=0,v=0))
outline=read_lines(url_outline)
url=sprintf("http://app-davinci.mars.asu.edu//Davinci/PSQLTESQueryTool?queryString=select %s from tes_4_jmars where ST_contains('%s',fp) and orbit_counter_keeper>%i and orbit_counter_keeper<%i and latitude between %f and %f and emission<%f and incidence<%f and tot_dust<%f and tot_ice<%f and lambert_albedo<%f and target_temperature>%f&displayHeader=1",query,outline[,1],ock_start,ock_end,minlat,maxlat,em_angle,in_angle,dust,waterice,albedo,temperature)
tes=load_csv(url)
################
# Extract information about emissivity, ock, ick, detector, tesx
N=length(tes.st_astext)
emi=create(N,143,1,format=float,start=0,step=0)
for(j=1;j<=N;j++) {
emi_tmp=strsplit(tes.emissivity[,j])
emi_tmp[,1]=delim(emi_tmp[,1],"{",1)
emi_tmp[,143]=delim(emi_tmp[,143],"}",1)
emi[j,,]=atof(emi_tmp)
}
ock=translate(tes.orbit_counter_keeper,x,y)
ick=translate(tes.instrument_time_count,x,y)
det=translate(tes.detector,x,y)
tesx=make_tesx()
tesx73=cat(tesx[,,9:35],tesx[,,65:110],z)
emi=translate(emi,y,z)
emi73=cat(emi[,,9:35],emi[,,65:110],z)
out=struct()
out.wkt=tes.st_astext
out.tesx=tesx
out.tesx73=tesx73
out.emi=emi
out.emi73=emi73
out.ock=ock
out.ick=ick
out.det=det
out.ock_ick_det=translate(cat(int(ock),int(ick),int(det),y),x,y)
ind=create(1,length(out.ock),1,start=1,step=1,format=int)
out.ock_ick_det=cat(ind,out.ock_ick_det,x)
out.minlat=minlat
out.maxlat=maxlat
if(minlat==-60 && maxlat==60) {
remove_struct(out,"minlat")
remove_struct(out,"maxlat")
}
# Display TES footprints overlapping THEMIS image
if(fexists(id+".cub")==1) {
geostruct=read_geo(id+".cub")
shape=struct()
shape.wkt=out.wkt
shape.string_feature=text(N)
shape.string_feature[,1:]="polygon"
geostruct.data=sstretch(geostruct.data[,,9])
overlap=shape2georaster(shape,geostruct,ignore=0,color=255//0//0,hsv=1)
display(overlap.data)
}
verbose=1
return(out)
}
define themis_atm_corr(minlib,atmlib,radius,power,iter,size){
if($ARGC==0) {
printf("\n")
printf("Use each overlying TES data to do THEMIS atmospheric correction automatically\n")
printf("Additionally, it will interpolate the atm component from TES data for gap regions\n")
printf("It will generate a fully atm corrected emissivity data for entire THEMIS image\n")
printf("$1: THEMIS id\n")
printf("$2: overlying TES data with necessary information from 'themis_tes_overlay'\n")
printf("Optional arguments:\n")
printf("minlib: Mineral library (default $DV_SCRIPT_FILES/themis_atm_support/THEMIS_ATM_mineral_library.hdf)\n")
printf("atmlib: 73 channel atmospheric library (default $DV_SCRIPT_FILES/themis_atm_support/TES_atmlib_73channel.hdf)\n")
printf("radius: radius for interpolation with weight_array\n")
printf("power: order for interpolation with weight_array\n")
printf("iterations: number of interations for interpolation running\n")
printf("size: square size for boxfilter to convolve atm component\n")
printf("Output: Corrected THEMIS IR cube\n")
printf("\n")
return(null)
}
verbose=0
id=$1
tes=$2
# mineral and atm spectral library
if(HasValue(minlib)==0) minlib=read($DV_SCRIPT_FILES+"/themis_atm_support/THEMIS_ATM_mineral_library.hdf")
if(HasValue(atmlib)==0) atmlib=read($DV_SCRIPT_FILES+"/themis_atm_support/TES_atmlib_73channel.hdf")
# radius and power designated for interpolated filling function
if(HasValue(radius)==0) radius=7
if(HasValue(power)==0) power=2
# iternations
if(HasValue(iter)==0) iter=0
# boxfilter size for atm component convolve
if(HasValue(size)==0) size=15
stime=syscall("date")[,1]
image=read_geo(id+".cub")
image.data[where min(image.data,axis=z)<-2]=-32768
geostruct=image
geostruct.data=sstretch(geostruct.data[,,9])
printf("\nTES surface emissivity calculation......\n")
spec=sma(tes.emi73,minlib,atmlib,wave1=200,wave2=1305,surface=1,nn=1)
surf_em=i2i(spec.rematm,from='tes73',to='themis')
shape=struct()
shape.wkt=text(1)
shape.string_feature=text(1)
shape.string_feature[,1]="polygon"
miny=1
maxy=image.h
if(HasValue(get_struct(tes,"maxlat"))!=0) {
miny=proj_geocoord(float(0.)//float(tes.maxlat),geostruct)
miny=int(round(miny[2]))
}
if(HasValue(get_struct(tes,"minlat"))!=0) {
maxy=proj_geocoord(float(0.)//float(tes.minlat),geostruct)
axy=int(round(maxy[2]))
}
N=dim(tes.emi73)[1]
atm_factor=create(image.w,maxy-miny+1,9,start=0,step=0,format=float)
mask_cont=create(image.w,maxy-miny+1,1,start=0,step=0,format=float)
for(i=1;i<=N;i++) {
shape.wkt[,1]=tes.wkt[,i]
mask=shape2georaster(shape,geostruct,ignore=0,color=1)
train_data=mask.data[,miny:maxy,1] * image.data[,miny:maxy,]
mask_cont=mask.data[,miny:maxy,1]+mask_cont
data=train_data[,,1:9]
data=thm.themis_emissivity(data,1//2//3//4//5//6//7//8//9,b1=3,b2=9)
data=data.emiss
ave_data=avg(data,axis=xy,ignore=0)
data [where (data<0.2)]=0
data [where (data>1.8)]=0
atm_trans=float(ave_data[,,1:9]/surf_em[i,,1:9])
atm_factor[where train_data!=0.]=atm_trans+atm_factor
}
printf("Start: %s\n\n", stime)
printf("Finish: %s\n\n", syscall("date")[,1])
atm_factor=atm_factor/mask_cont
emi_data=thm.themis_emissivity(image.data[,,1:9],1//2//3//4//5//6//7//8//9,b1=3,b2=9)
emi_data=emi_data.emiss
emi_data [where (emi_data<0.2)]=0
emi_data [where (emi_data>1.8)]=0
output=struct()
output.atm_opacity=-ln(atm_factor)
output.surf_emi = emi_data / atm_factor
# filling the gap between TES pixels
wt=weight_array(radius,power)
for (j=1;j<=9;j++) {
atm_factor[,,j]=kfill(atm_factor[,,j],wt,iterations=iter)
}
atm_factor=boxfilter(atm_factor,size=size,ignore=0)
output.atm_opacity_fill=-ln(atm_factor)
boundary=image.data[,,9]
boundary[where min(boundary,axis=z)!=-32768]=1
boudnary[where min(boundary,axis=z)==-32768]=0
atm_factor=atm_factor * boundary
output.surf_emi_fill = emi_data / atm_factor
write(output,id+'_autorematm.hdf',hdf,force=1)
out=image
out.data=output.surf_emi_fill
write_geotiff(out,id+'_autorematm_gapfill.tiff',force=1)
return(output)
}
define tes_fp(){
if($ARGC==0){
printf("\Get coordinates and topography information of TES data overlapping THEMIS image\n")
printf("\$1: THEMIS ID\n")
printf("\$2: Result from themis_tes_overla() function\n")
printf("\output: Structure of all TES data with all necessary information\n\n")
return(null)
}
id=$1
tes=$2
geostruct=read_geo(id+".cub")
# Find latitude and longitude of TES polygon
N=length(tes.wkt)
clon=create(1,1,1,start=0,step=0,format=float)
clat=create(1,1,1,start=0,step=0,format=float)
ullon=create(1,1,1,start=0,step=0,format=float)
ullat=create(1,1,1,start=0,step=0,format=float)
lrlon=create(1,1,1,start=0,step=0,format=float)
lrlat=create(1,1,1,start=0,step=0,format=float)
for(i=1;i<=N;i++) {
tes_poly=strsplit(tes.wkt[,i],",")
tes_poly[,1]=delim(tes_poly[,1],"((",2)
V=length(tes_poly)
tes_lon=delim(tes_poly[,1:V-1]," ",1)
tes_lat=delim(tes_poly[,1:V-1]," ",2)
tes_clon=float(avg(atof(tes_lon)))
tes_clat=float(avg(atof(tes_lat)))
tes_ullon=min(atof(tes_lon))
tes_ullat=max(atof(tes_lat))
tes_lrlon=max(atof(tes_lon))
tes_lrlat=min(atof(tes_lat))
clon=cat(clon,tes_clon,y)
clat=cat(clat,tes_clat,y)
ullon=cat(ullon,tes_ullon,y)
ullat=cat(ullat,tes_ullat,y)
lrlon=cat(lrlon,tes_lrlon,y)
lrlat=cat(lrlat,tes_lrlat,y)
}
clon=clon[,2:,]
clat=clat[,2:,]
ullon=ullon[,2:,]
ullat=ullat[,2:,]
lrlon=lrlon[,2:,]
lrlat=lrlat[,2:,]
# Transfer lat/lon to x/y (using array instead of single number)
tes_cxy=proj_geocoord(clon//clat,geostruct)
tes_cxy=int(round(tes_cxy))
tes_ulxy=proj_geocoord(ullon//ullat,geostruct)
tes_ulxy=int(round(tes_ulxy))
tes_lrxy=proj_geocoord(lrlon//lrlat,geostruct)
tes_lrxy=int(round(tes_lrxy))
# Topography information
topo=get_map(geostruct,map="mola_hrsc_blend_topo_v2",trim=1)
out=struct()
out.wkt=tes.wkt
out.center=tes_cxy
out.ul=tes_ulxy
out.lr=tes_lrxy
out.topo=topo.data
out.tesx=tes.tesx
out.tesx73=tes.tesx73
out.emi=tes.emi
out.emi73=tes.emi73
out.ock=tes.ock
out.ick=tes.ick
out.det=tes.det
out.ock_ick_det=tes.ock_ick_det
# Display TES footprints overlapping THEMIS image
shape=struct()
shape.wkt=out.wkt
shape.string_feature=text(N)
shape.string_feature[,1:]="polygon"
geostruct.data=sstretch(geostruct.data[,,9])
overlap=shape2georaster(shape,geostruct,ignore=0,color=255//0//0,hsv=1)
display(overlap.data)
return(out)
}
define refine_tes_fp(dis,elv){
if($ARGC==0){
printf("\Find TES data that meeting the geography and elevation constrains used for training area\n")
printf("\$1: THEMIS ID\n")
printf("\$2: Results from tes_fp with coordinate information\n")
printf("\$3: x1//y1, upper left corner of THEMIS targeted study area\n")
printf("\$4: x2//y2, lower right corner of THEMIS targeted study area\n")
printf("\$5: x1//y1, upper left corner of training area\n")
printf("\$6: x2//y2, lower right corner of training area\n")
printf("\dis: optional, TES footpoint and THEMIS ROI distance in y\n")
printf("\elv: optional, TES footpoint and THEMIS ROI elevation difference\n")
printf("\output: a TES geostruct file meeting the constrains\n\n")
return(null)
}
id=$1
tes=$2
roi_ul=$3
roi_lr=$4
training_ul=$5
training_lr=$6
if(HasValue(dis)==0) dis=1000
if(HasValue(elv)==0) elv=1000
geostruct=read_geo(id+".cub")
geostruct.data=sstretch(geostruct.data[,,9])
# Get targeted study area lat/lon coordinate
roi_ullatlon=proj_geocoord(roi_ul,geostruct,inverse=0)
roi_lrlatlon=proj_geocoord(roi_lr,geostruct,inverse=0)
roi_ullon=roi_ullatlon[1]
roi_ullat=roi_ullatlon[2]
roi_lrlon=roi_lrlatlon[1]
roi_lrlat=roi_lrlatlon[2]
if(roi_ullon<0) roi_ullon=360.0+roi_ullon
if(roi_lrlon<0) roi_lrlon=360.0+roi_lrlon
roi_polygon=sprintf("POLYGON((%f %f,%f %f,%f %f,%f %f,%f %f))",roi_ullon,roi_ullat,roi_lrlon,roi_ullat,roi_lrlon,roi_lrlat,roi_ullon,roi_lrlat,roi_ullon,roi_ullat)
roi=struct()
roi.wkt=text(1)
roi.wkt[,1]=roi_polygon
roi.string_feature=text(1)
roi.string_feature[,1]="polygon"
# overlap=shape2georaster(roi,geostruct,ignore=0,color=0//255//0,hsv=1)
# display(overlap.data)
roi_cx=(roi_ul[1]+roi_lr[1])/2
roi_cy=(roi_ul[2]+roi_lr[2])/2
roi_elv=avg(tes.topo[roi_ul[1]:roi_lr[1],roi_ul[2]:roi_lr[2],])
#roi_elv=tes.topo[roi_cx,roi_cy,]
# Find TES data meeting the geography and elevation criteria from training area
N=length(tes.wkt)
wkt=""
emi=create(1,1,143,start=0,step=0,format=float)
emi73=create(1,1,73,start=0,step=0,format=float)
center=create(2,1,1,start=0,step=0,format=int)
ul=create(2,1,1,start=0,step=0,format=int)
lr=create(2,1,1,start=0,step=0,format=int)
ock=create(1,1,1,start=0,step=0,format=int)
ick=create(1,1,1,start=0,step=0,format=int)
det=create(1,1,1,start=0,step=0,format=int)
for(i=1;i<=N;i++) {
diff_dis=abs(tes.center[2,i]-roi_cy)
#diff_x=tes.center[1,i]-roi_cx
#diff_y=tes.center[2,i]-roi_cy
#diff_dis=hyp(diff_x,diff_y)
tes_elv=avg(tes.topo[tes.ul[1,i]:tes.lr[1,i],tes.ul[2,i]:tes.lr[2,i],])
#tes_elv=tes.topo[tes.center[1,i],tes.center[2,i],]
diff_elv=abs(roi_elv-tes_elv)
if(diff_dis
if(diff_elv
if((tes.center[1,i]>=training_ul[1]) && (tes.center[2,i]>=training_ul[2]) && (tes.center[1,i]<=training_lr[1]) && (tes.center[2,i]<=training_lr[2])) {
wkt=cat(wkt,tes.wkt[,i],y)
center=cat(center,tes.center[,i],y)
ul=cat(ul,tes.ul[,i],y)
lr=cat(lr,tes.lr[,i],y)
emi=cat(emi,tes.emi[i],x)
emi73=cat(emi73,tes.emi73[i],x)
ock=cat(ock,int(tes.ock[i]),x)
ick=cat(ick,int(tes.ick[i]),x)
det=cat(det,int(tes.det[i]),x)
}
}
}
}
num=length(wkt)
if(num<=1) {
printf("\nNo TES data found meeting the search criteria!\n")
return(null)
}
new_wkt=wkt[,2:]
new_center=center[,2:]
new_ul=ul[,2:]
new_lr=lr[,2:]
new_emi=emi[2:]
new_emi73=emi73[2:]
new_ock=ock[2:]
new_ick=ick[2:]
new_det=det[2:]
out=struct()
out.wkt=new_wkt
out.center=new_center
out.ul=new_ul
out.lr=new_lr
out.emi=new_emi
out.emi73=new_emi73
out.ock=new_ock
out.ick=new_ick
out.det=new_det
out.tesx=tes.tesx
out.tesx73=tes.tesx73
out.topo=tes.topo
out.tes_label=translate(cat(new_ock,new_ick,new_det,y),x,y)
shape=struct()
shape.wkt=out.wkt
shape.string_feature=text(length(out.wkt))
shape.string_feature[,1:]="polygon"
overlap=shape2georaster(shape,geostruct,ignore=0,color=255//0//0,hsv=1)
display(overlap.data)
out.roi=struct()
out.roi.ullatlon=roi_ullatlon
out.roi.lrlatlon=roi_lrlatlon
out.roi.ulxy=roi_ul
out.roi.lrxy=roi_lr
return(out)
}
define themis_rematm(minlib,atmlib,tes1,tes2){
if($ARGC==0) {
printf("\n")
printf("Atmospheric correction of THEMIS targeted study area\n")
printf("$1: THEMIS id\n")
printf("$2: TES footprints used as training area\n")
printf("minlib: Mineral library (default $DV_SCRIPT_FILES/themis_atm_support/THEMIS_ATM_mineral_library.hdf)\n")
printf("atmlib: 73 channel atmospheric library (default $DV_SCRIPT_FILES/themis_atm_support/TES_atmlib_73channel.hdf)\n")
printf("optional:\n")
printf("tes1=start orbit index for tes\n")
printf("tes2=end orbit index for tes\n")
printf("Output: Corrected THEMIS data\n")
printf("\n")
return(null)
}
id=$1
tes=$2
if(HasValue(minlib)==0) minlib=read($DV_SCRIPT_FILES+"/themis_atm_support/THEMIS_ATM_mineral_library.hdf")
if(HasValue(atmlib)==0) atmlib=read($DV_SCRIPT_FILES+"/themis_atm_support/TES_atmlib_73channel.hdf")
if(HasValue(tes1)==0) tes1=1
if(HasValue(tes2)==0) tes2=length(tes.ock)
image=read_geo(id+".cub")
geostruct=image
geostruct.data=sstretch(geostruct.data[,,9])
tes_avg=avg(tes.emi73[tes1:tes2],axis=x)
spec=sma(tes_avg,minlib,atmlib,wave1=200,wave2=1305,surface=1,nn=1)
surf_em=i2i(spec.rematm,from='tes73',to='themis')
N=dim(tes.emi73)[1]
shape=struct()
shape.wkt=tes.wkt
shape.string_feature=text(N)
shape.string_feature[,1:]="polygon"
roi_ul=tes.roi.ulxy
roi_lr=tes.roi.lrxy
mask=shape2georaster(shape,geostruct,ignore=0,color=1)
training_area=mask.data[,,1] * image.data
tes_em=surf_em[,,1:9]
ave_data=training_area[,,1:9]
ave_data=thm.themis_emissivity(ave_data, 1//2//3//4//5//6//7//8//9,b1=3,b2=9)
ave_data=avg(ave_data.emiss,axis=xy,ignore=0.0)
data=image.data[,,1:9]
data=thm.themis_emissivity(data, 1//2//3//4//5//6//7//8//9,b1=3,b2=9)
data=data.emiss
data [where (data<0.2)]=0
data [where (data>1.8)]=0
atm_trans=float(ave_data[,,1:9]/tes_em)
thmx=make_band(themis=1,cm=1,descend=1)
out=struct()
out.surf=data/atm_trans
out.atm=-ln(atm_trans)
out.xaxis=thmx[,,1:9]
write(out,id+'_semirmatm.hdf',hdf,force=1)
return(out)
}