Monday, November 30, 2015

Standardized Precipitation-Evapotranspiration Index

The SPEI fulfils the requirements of a drought index since its multi-scalar character enables it to be used by different scientific disciplines to detect, monitor and analyze droughts. Like the sc-PDSI and the SPI, the SPEI can measure drought severity according to its intensity and duration, and can identify the onset and end of drought episodes. The SPEI allows comparison of drought severity through time and space, since it can be calculated over a wide range of climates, as can the SPI. Moreover, Keyantash and Dracup (2002) indicated that drought indices must be statistically robust and easily calculated, and have a clear and comprehensible calculation procedure. All these requirements are met by the SPEI. However, a crucial advantage of the SPEI over other widely used drought indices that consider the effect of PET on drought severity is that its multi-scalar characteristics enable identification of different drought types and impacts in the context of global warming.

SPEI详细的介绍请见Vicente-Serrano SM, et al., 2010。随论文，Vicente-Serrano开发了SPEI Calculator程序，其中PET应用仅包含温度参数的计算方法(Thornthwaite, 1948)，这个程序自由下载，使用说明参考spei_manual_en.pdf。

应对SPEI批量操作需求，我在不掌握SPEI计算过程的基础上，以Matlab调用SPEI Calculator程序实现批量SPEI指数计算，具体情况见附件中SPEI_Cal.m和SPEI_Batch.m，代码运行流畅，调用返回结果正确。

图 1~3分别表示1-、3-、6-、9-、12-month尺度SPEI指数，Orig-12-month SPEI是测试数据附带的测试结果。

图 1

图 2

图 3

如图 3所示，上下两图干湿交替情况基本相同。不过，仔细查看两组数据，我们发现数据之间仍存在微小的差异，比如1901年1月12-month的SPEI在上下两图分别为1.152和1.164104，原因可能在于SPEI Calculator运行环境的差异，但不表示本地运行的前者是错误结果。

SPEI干旱/湿润等级划分(Yu, et al., 2014)


Categories	SPEI values

Extremely dryness	Less than - 2
Severe dryness	-1.99 to -1.5
Moderate dryness	-1.49 to -1.0
Near normal	-1.0 to 1.0
Moderate wetness	1.0 to 1.49
Severe wetness	1.50 to 1.99
Extremely wetness	More than -2

附上测试代码及数据（SPEI.rar）。

附上栅格处理参考代码（SPEIcal.rar）。

References

[1] SPEI.

[2] SPEI Calculator.

[3] Vicente-SerranoSM, Beguería S, López-Moreno JI (2010) A multiscalar drought index sensitive to global warming: the standardized precipitation evapotranspiration index. Journal of Climate, 23, 1696-1718.

[4] Thornthwaite CW (1948) An approach toward a rational classification of climate. Geographical Review, 55-94.

[5] How to run .exe file from matlab.

[6] Yu M, Li Q, Hayes MJ, Svoboda MD, Heim RR (2014) Are droughts becoming more frequent or severe in China based on the standardized precipitation evapotranspiration index: 1951–2010? International Journal of Climatology, 34, 545-558.

Matlab: Rendering Automatically

Summary

全自动二维图像渲染程序，将单波段图像按照唯一数值（Unique Value）和区间数值（Classified Value）进行RGB渲染。

Settings

文件Info.xlsx记录渲染配置的颜色，分别是唯一数值（Fig. 1）和区间数值（Fig. 2），第一列标示代码的含义（中英文皆可，服务于理解），不做运算处理，之后的2~3或2~4列参与运算。在区间数值运算中，区间遵循（第二列，第三列】规则对应第四列RGB。

Fig. 1

Fig. 2

Fig. 3是代码运行结果，Fig. 4是手动渲染结果，两者别无二致。

Fig. 3

Fig. 4

注意：当代码渲染结果放入ArcGIS中，它会被自动拉伸，这并不是我们想要的情况，只要在属性（Properties）中设置Stretch Type: None，并取消勾选Apply Gamma Stretch，见Fig. 5。

Fig. 5

Rendering.rar.

Rending and Exporting to Png format images.

Friday, November 27, 2015

Math: Tips

[1] 绝对变化与相对变化？

Wednesday, November 25, 2015

Matlab: Quadratic Fit

周期性变量的趋势不能用线性回归进行拟合，周期性波动不符合线性特征，可以考虑应用一元二次曲线拟合，例如月降水量在一年内的变化。

图 1是某地多年1~12月平均降水量的变化情况，一年12个月正好是月降水量的完整变化周期，尝试一元二次曲线拟合这个周期性的变化趋势，以R²表示方差解释率，p-value表示回归模型的显著性。

Fig. 1

图 2显示R²为0.5602，表示约有56%的方差可以被解释，比例还比较高。p-value为0.0248，若以α=0.05为显著水平，则p-value＜α，通过显著性检验，该模型成立。

Fig. 2

Mat2Polyfit.rar.

栅格运算（单波段、多波段）.

Tuesday, November 24, 2015

Matlab: 不同Size图片镶嵌

开贴介绍不同Size图片的镶嵌方法，不同Size的图片是指图片的行列不相同的多个图片。

以两个不同Size的图片为例，运用代码实现图片在右侧或下方自动调整后的镶嵌操作，1.jpg有817行815列，2.jpg有799行800列。我们的目标是实现图 1的两种效果。镶嵌过程中1.jpg的Size不变，但2.jpg依情况有对应的调整。

图 1

函数imresize()操作2.jpg，依据镶嵌的行或列调整2.jpg。以水平镶嵌（keyword控制）为例，根据1.jpg有817行而得调整之后的2.jpg的伸缩Scale是1.0225，伸缩后Size是817×819。代码处理结果如图 2和3，毫无违和。

图 2

图 3

MatMerge.rar.

MatMerge_2.rar.上下左右批量镶嵌.2018-03-02.

Tool: 网易邮箱插入图像签名的方法

最终方法是：

上传图片至网易相册（或者其他更容易获取的地方，如公开的博客），获取该图片地址。
打开网易邮箱，点写信，在写信正文栏插入该图片。
点击选取图片，调整至中意的签名用大小尺寸，然后复制该图片。
点右上角“设置”，然后选“签名设置”；在签名设置页面下点“添加签名”。
在添加签名的正文部分点右键粘贴图片，即能把该图片放入签名档。
加入签名档名称及签名的其他内容，保存即可。

Matlab: Enhanced Mosaic Operation

Summary

若干图片还是要求具有同样的尺寸，并含有相同的波段数量。

Code Specification

代码前端设定这些内容:

N=2;定义输出列数
color=[255, 255, 255];定义输出背景颜色
blkimg.width=0.04;定义插入图片之间水平方向间距，归一化为图片的width
blkimg.height=0.04;定义插入图片之间垂直方向间距，归一化为图片的height
edgimg.width=0;定义图片与边缘之间水平方向间距，归一化为图片的width
edgimg.height=0;定义图片与边缘之间垂直方向间距，归一化为图片的height

如果blkimg.width、blkimg.height、edgimg.width、edgimg.height四变量全部赋值为0，此时结果将与下面代码结果一致，不留出任何空白，见Fig. 1。

% Created by LI Xu
% Version 1.0
% 25 October, 2014

% Modified by LI Xu
% Version 1.1
% November 24, 2015
% Add Source and destination directory

% Mosaicing Images with the same
% rows, columns and banks.
% No blank edges.

clear;
clc;

% Set the number of columns for Output M*N
N=2;

% Source Directory
SouDir='./input';
% Destination Directory
DesDir='./output';

% Images
files=dir([SouDir, '/*.png']);
num_files=length(files);
M=ceil(num_files/N);

% Image Info
img=imread([SouDir, '/', files(1).name]);
rows=size(img, 1);
cols=size(img, 2);
bands=size(img, 3);

cellstack=cell(M, 1);
N_flag=1;
temp_img=zeros(rows, cols*N, bands)+255;
temp_img=uint8(temp_img);
for ii=1:num_files
    imgfile=files(ii).name;
    in_img=imread([SouDir, '/', imgfile]);
    
    % Position
    if N_flag > N
        N_flag=1;
        temp_img=zeros(rows, cols*N, bands)+255;
        temp_img=uint8(temp_img);
    end
    cellrow=ceil(ii./N);
    row_position=1:rows;
    col_position=1+cols*(N_flag-1):cols*N_flag;
    temp_img(row_position, col_position, :)=in_img;
    %imshow(temp_img);
    cellstack{cellrow, 1}=temp_img;
    
    N_flag=N_flag+1;
end


% Output Image
otImage=zeros(rows*M, cols*N, bands);
for row=1:M
    row_range=1+rows*(row-1):rows*row;
    otImage(row_range, :, :)=cellstack{row, 1};
end

otImage=uint8(otImage);
imshow(otImage);
pause(2);

otname=['NoEdge_R', num2str(M), 'C', num2str(N), '.png'];
otpath=[DesDir, '/', otname];
imwrite(otImage, otpath);
close all
disp('END');

Fig. 1

其他时候可以是这样的：

Fig. 2

Fig. 3

MatMosaic.rar.

MatMosaic1.2.rar.灰度与彩色图片混搭.

不同尺寸的图片.

References

[1] Matlab: 图片镶嵌.

Monday, November 23, 2015

Matlab: Insert Text

Summary

Matlab如何在图片上更灵活地插入文字是一个不断发展的主题，立足于自身需求及对代码的理解，这一问题今天有一新的回答。当然，目前Matlab在图片中插入文字仅仅支持ASCII编码，也就是不支持中文字符。

Code Specification

代码重点在设置待插入文字的属性。举例：

stratts_1.fontname='Arial';
stratts_1.fontsize=0.08; 文字大小，归一化为图片height的比例
stratts_1.fontcolor=[255, 215, 0]; 文字颜色
stratts_1.xscale=0.05; X轴起始位置，归一化为图片width的比例
stratts_1.yscale=0.05; Y轴起始位置，归一化为图片height的比例
stratts_1.rowinterval=0.2; 多行文字的行间距，归一化为图片height的比例
stratts_1.Opacity=1; 文字透明度

此处，代码从.xslx文件中录入待插入文字，属性fontsize、xscale、yscale、rowinterval若大于1则认为是绝对数值，增加代码的可移植性。以下Fig.1~4是一些插入文字的结果：

Fig. 1

Fig. 2

Fig. 3

Fig. 4

这个小程序还有不足之处，未来继续改进。

MatInserText.rar.

[2023-10-16 Update]Insert Text which was from the file name.

References

[1] 文字加入图片.

Sunday, November 22, 2015

Matlab: Pointy Ends for Colorbars

Cbarrow

The cbarrow function places triangle-shaped endmembers on colorbars to indicate that data values exist beyond the extents of the values shown in the colorbar.

This function works by creating a set of axes atop the current figure and placing patch objects on the new axes. Thus, editing a figure after calling cbarrow may cause some glitches. Therefore, it is recommended to call cbarrow last when creating plots.

This function only works once per figure. If you have multiple subplots, you can only use it once, and you'll have to call cbarrow last. Also, editing plots after calling cbarrow can sometimes be a bit glitchy.

Fig. 1

Fig. 2

Cbarrow函数要求Matlab版本至少是2014b才能运行。

Cabrf

Fig. 3

Attachment.

References

[1] cbarrow: pointy ends for colorbars - home page.

[2] ColorBrewer: Attractive and Distinctive Colormaps.

[3] cbarf - home page.

Python: 中英文字符插入图片

Summary

中文字符批量插入图片，中文属性配置在Info.xlsx文件，包括所在图片名称关键字、插入字符、字体大小、字体名称、X相对位置以及Y相对位置。

Note

中文字符编码需要格外注意（冷僻汉字不一定可以顺利插入），在代码中要有编码声明，这里：#-*- coding:utf8 -*-，知会计算机中文字符的编码格式。
Xls文件一定是utf-8编码格式，该文件编码调整方法见这里。
文本PyAddText.py本身也必是utf-8编码格式，该文件编码调整方法见这里。
字体名称一定对应着字体文件名称，如字体Arial对应字体文件arial.ttf，在属性字体中设置‘arial’，如平常Office显示之Arial则出现错误。

满足以上条件之后，代码可以正确插入中文。不过，每次运行都会发出警告信息，这个留在以后再解决。

"Warning (from warnings module): File "C:\Python27\ArcGISx6410.2\lib\site-packages\PIL\ImageFont.py", line 273 if ext and walkfilename == ttf_filename: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal "
-from SHELL Window

Fig. 1

Fig. 2

Fig. 3

#-*- coding:utf8 -*-
# Created by LI Xu
# Version 1.0
# November 22, 2015

# Description:
# Insert Chinese Characters on pictures

# If you have any question about this code,
# please do not hesitate to contact me via E-mail: 
# jeremy456@163.com

# Blog:
# http://blog.sciencenet.cn/u/lixujeremy
# http://lixuworld.blogspot.com/


import sys
import math
import csv
import os

from PIL import Image,ImageDraw,ImageFont
import datetime
import time
begintime=time.strftime("%Y-%m-%d %H:%M:%S")
print 'Time Begin:'+begintime
starttime = datetime.datetime.now()

#~#------------------------------------------------------------------
def GetXlsRows(xlspath, indexes, records, offset):
    import xlrd
    workbook=xlrd.open_workbook(xlspath)
    # loop for every single sheet
    import math
    for sheetno in indexes:
        #print sheetno
        try:
            worksheet=workbook.sheet_by_index(sheetno)
            records=GetRecords(worksheet, offset, records)
        except:
            continue
    return records
#~#------------------------------------------------------------------
def GetRecords(worksheet, offset, records):
    for i, row in enumerate(range(worksheet.nrows)):
        if i <offset:
            continue
        r=[]    
        for j, col in enumerate(range(worksheet.ncols)):
            r.append(worksheet.cell_value(i, j))
        records.append(r)
    return records
#~#------------------------------------------------------------------
def GetFileSubStr(SouDir, flagstr):
    import os, glob
    allfiles=os.listdir(SouDir)
    filepath=[]
    for file in allfiles:
        if flagstr in file:
            filepath=os.path.join(SouDir, file)
            return filepath

#~ #----------------------------------------------------------------------  

#~ #----------------------------------------------------------------------
def InsertTextOnImg(SouDir, DesDir, record):
    flagstr=record[0]
    filepath=GetFileSubStr(SouDir, flagstr)
    if not filepath:
        print 'No such file: '+flagstr
        return
    InWord=record[1]
    FontSize=int(record[2])
    FontName=record[3]+'.ttf'
    Scale_1=record[4]
    Scale_2=record[5]
    #Output filepath
    filename=filepath.split('\\')
    filename=filename[len(filename)-1]
    otpath=os.path.join(DesDir, filename)
    print otpath
    imgfile=Image.open(filepath)
    imgSize=imgfile.size
    draw=ImageDraw.Draw(imgfile)
    font=ImageFont.truetype(FontName, FontSize)
    draw.text(((imgSize[0])*Scale_1, (imgSize[1])*Scale_2), InWord, (0,0,0), font=font)
    imgfile.save(otpath, 'PNG')

#~ #---------------------------------------------------------------------- 

#MAIN EXTRANCE
#Get the current working directory
WDir=os.getcwd()
#Source Directory
SouDir=os.path.join(WDir, 'input')
#Destination Directory
DesDir=os.path.join(WDir, 'output')
#Info
infopath='Info.xlsx'
#Change this depanding on how many header rows are present
#set to 0 if you want to include the header data.
offset=1
records=[]
index=[0, 1, 2]
info=GetXlsRows(infopath, index, records, offset)

#Loop for every single record
for ii in range(0, len(info)):
    rocd=info[ii]
    print ii+1,
    #print rocd
    InsertTextOnImg(SouDir, DesDir, rocd)


print 'Time Begin:'+begintime
starttime = datetime.datetime.now()
print "END"
endtime=time.strftime("%Y-%m-%d %H:%M:%S")
print 'Time End:'+endtime
finishtime=datetime.datetime.now()
timespan=(finishtime-starttime).seconds
timespan='%f' %timespan
print 'Time Span:'+timespan+' s'

print '************************************'

Attachment.

References

[1] Python+Pillow: 中文字符插入图片.

Saturday, November 21, 2015

Tool: MS Office

[1] How to Encode an Excel File to UTF-8 or UTF-16?

[2] [Endnote]Chinese Standard GB/T7714 (numeric).

[3] 文章里“摄氏度”符号的输入方法.

[4] Insert a nonbreaking minus sign? Control+Shift+-.

[5] Cannot use EZ-Activator because your licensing details are corrupt. [Office 2010 Toolkit]

[6] 外国人名字中间的点怎么打.

[7] Word2010点击保存按钮无法保存，另存为按钮跳不出对话框.

[8] 文字底纹.

[9] Win10环境下设置Office背景豆沙绿【前两步即可】.

Friday, November 20, 2015

Matlab: Converison from Arc/Info Grid to GeoTiff Format

Summary

笔记：记录ANUSPLIN LAPGRD输出.grd文件转换为通用文件Geotiff之方法。

Note

The codes in this case is to be oriented to the files associated with Lat/Long coordinate system only. For others, I would code in the future probably.

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

clear;
clc;

% Source Directory
SouDir='./Input';
% Destination Directory
DesDir='./Output';

% files
files=dir([SouDir, '/*.grd']);

for ii=1:length(files)
    filename=files(ii).name;
    filepath=[SouDir, '/', filename];
    
    % Output file fullpath
    outname=[filename(1:end-3), 'tif'];
    outpath=[DesDir, '/', outname];
    
    % Conversion
    Gonvert2GTiff(filepath, outpath);
    
    disp([num2str(ii), ':', outpath]);
end


disp('***************************************');

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

function Gonvert2GTiff(grdpath, gtiffpath);

    % Get the properties of this .grd file
    GeoInfo=GetProperty(grdpath);
    
    % Get the values of this .grd file;
    MatVal=GetValues(grdpath);
    MatVal=reshape(MatVal, GeoInfo.NCOLS, GeoInfo.NROWS);
    MatVal=MatVal';
    
    
    % Spatial extent
    stLong=GeoInfo.XLLCORNER;
    edLong=stLong+GeoInfo.CELLSIZE*GeoInfo.NCOLS;
    stLat=GeoInfo.YLLCORNER;
    edLat=stLat+GeoInfo.CELLSIZE*GeoInfo.NROWS;
    
    % Write a Data Referenced to Geographic Coordinates
    R=georasterref('RasterSize', [GeoInfo.NROWS, GeoInfo.NCOLS], 'Latlim', [stLat, edLat], 'Lonlim', [stLong, edLong], 'ColumnsStartFrom', 'north');
    
    geotiffwrite(gtiffpath, MatVal, R);



end

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Get values of a .grd file

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

function Output=GetValues(grdpath)
    
    % Read values
    fid=fopen(grdpath);
    
    tline=fgets(fid);
    count=1;
    Output=[];
    while ischar(tline)
        
        if count>6
            Vals=str2num(tline);
            Output=[Output; Vals'];
        end       
        tline=fgets(fid);
        count=count+1;
    end
    
     fclose(fid);
    fclose('all');


end

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy


function GeoInfo=GetProperty(grdpath)
    

    fid=fopen(grdpath);
    tline=fgets(fid);
    for ii=1:6
        [property, value]=strread(tline, '%s %f');
        str_exe=['GeoInfo.', upper(property{1}), '=', num2str(value), ';'];
        eval(str_exe);
        tline=fgets(fid);
    end
    
    fclose(fid);
    fclose('all');
    
end

References

[1] 实现气象数据的空间分布.

Matlab: Masking the Imagery

Matlab操作非常简单，mask文件的背景也是目标影像的背景，以mask影像的背景索引为目标影像的索引赋值为背景值，结果保存在新文件。

待掩膜的影像数据应与Mask.tif具有一致的行列及坐标、空间分辨率信息，如Fig. 1。可以在代码中配置输出影像数据的背景数值，BackValue=-99;。

Fig. 1

% Created by LI Xu
% Version 1.0
% 30 September, 2015

% Description:
% Masking the imageries.

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

clear;
clc;


% Mask
maskpath='mask.tif';
[mask, geo]=geotiffread(maskpath);
info=geotiffinfo(maskpath);
index_invalid=find(mask==0);

% Sourse Directory
SouDir='./new';
% Destination Directory
DesDir='./snow';

% All files to be masked
files=dir([SouDir, '/*.tif']);
% BackGround Value
BackValue=-99;

for ii=1:length(files)
    filename=files(ii).name;
    inpath=[SouDir, '/', filename];
    
    image=imread(inpath);
    image(index_invalid)=BackValue;
    
    % Output file path
    otpath=[DesDir, '/', filename];
    geotiffwrite(otpath, image, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp([num2str(ii), ':', otpath]);
    
end

disp('********************************');

IDL: 删除影像中的波段

问题：删除影像中的若干波段，并生成新影像文件？

方法：IDL。练习数据中，test.tif包含13个波段，删除其中的1~N（小于13）个波段，得到结果result.tif。

;+
; :Author: LiXu
;-


;;Created by LI Xu
;;Version 1.0
;;19 October, 2014

pro Delete_Bands_Batch

  print, 'Time Begin:'+' '+systime()
  begintime=SYSTIME(1)
  
  ;;Source Directory
  SouDir='E:\Tools\rad\Rad'
  ;;Destination Directory
  DesDir='E:\Tools\rad\NEW'
  ;;The Bands to Delete
  band=[1]
  
  
  ;;Retrieve all files
  filespath=file_search(SouDir, '*.tif', count = num_files, /test_regular)
  for ii=0,  num_files-1 do begin
    imgpath=filespath(ii)
    ;;otImagepath
    filename=strsplit(imgpath, '\', count=str_count, /Extract)
    filename=filename(str_count-1)
    otpath=DesDir+'\'+filename
    
    Delete_Bands, imgpath, otpath, band
    
  endfor
  
  
  

  print, 'End!'
  print, 'Time End:'+' '+systime()
  endtime=systime(1)
  timespan=endtime-begintime
  print, 'Time Span:'+' '+string(timespan)+' s'







end

;+
; :Author: LiXu
;-

;;Created by LI Xu
;;Version 1.0
;;18 October, 2014

;;Description:
;;Delete several bands from the input Image,
;;and generate a new Image.

pro Delete_Bands, infile, otfile, numoffile

  ;;Debug Block
  ;;infile='F:\Tools\Practice\test.tif'
  ;;otfile='F:\Tools\Practice\result.tif'
  ;;numoffile=[1]
  
  ;;Code Body
  imgfile=read_tiff(infile, GeoTIFF=GEOVAR)
  filesize=size(imgfile)
  Num_Bands=filesize(1)
  Bands=indgen(13)
  
  cols=filesize(2)
  rows=filesize(3)
  
  ;;The Bands to Delete
  Num_Dlet=N_Elements(numoffile)
  ;;The Bands to Save
  Num_Save=Num_Bands-Num_Dlet
  Remove, numoffile-1, Bands
  
  newfile=make_array(Num_Save, cols, rows, /float)
  for ii=0, Num_Save-1 do begin
    newfile[ii, *, *]=imgfile[Bands[ii], *, *]

  endfor

  ;;Write
  write_tiff, otfile, newfile, GEOTIFF=GEOVAR, /float





end

pro remove,index, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v14, $
  v15, v16, v17, v18, v19, v20, v21, v22, v23, v24, v25
  ;+
  ; NAME:
  ;       REMOVE
  ; PURPOSE:
  ;       Contract a vector or up to 25 vectors by removing specified elements
  ; CALLING SEQUENCE:
  ;       REMOVE, index, v1,[ v2, v3, v4, v5, v6, ... v25]
  ; INPUTS:
  ;       INDEX - scalar or vector giving the index number of elements to
  ;               be removed from vectors.  Duplicate entries in index are
  ;               ignored.    An error will occur if one attempts to remove
  ;               all the elements of a vector.     REMOVE will return quietly
  ;               (no error message) if index is !NULL or undefined.
  ;
  ; INPUT-OUTPUT:
  ;       v1 - Vector or array.  Elements specifed by INDEX will be
  ;               removed from v1.  Upon return v1 will contain
  ;               N fewer elements, where N is the number of distinct values in
  ;               INDEX.
  ;
  ; OPTIONAL INPUT-OUTPUTS:
  ;       v2,v3,...v25 - additional vectors containing
  ;               the same number of elements as v1.  These will be
  ;               contracted in the same manner as v1.
  ;
  ; EXAMPLES:
  ;       (1) If INDEX = [2,4,6,4] and V = [1,3,4,3,2,5,7,3] then after the call
  ;
  ;               IDL> remove,index,v
  ;
  ;       V will contain the values [1,3,3,5,3]
  ;
  ;       (2) Suppose one has a wavelength vector W, and three associated flux
  ;       vectors F1, F2, and F3.    Remove all points where a quality vector,
  ;       EPS is negative
  ;
  ;               IDL> bad = where( EPS LT 0, Nbad)
  ;               IDL> if Nbad GT 0 then remove, bad, w, f1, f2, f3
  ;
  ; METHOD:
  ;       If more than one element is to be removed, then HISTOGRAM is used
  ;       to generate a 'keep' subscripting vector.    To minimize the length of
  ;       the subscripting vector, it is only computed between the minimum and
  ;       maximum values of the index.   Therefore, the slowest case of REMOVE
  ;       is when both the first and last element are removed.
  ;
  ; REVISION HISTORY:
  ;       Written W. Landsman        ST Systems Co.       April 28, 1988
  ;       Cleaned up code          W. Landsman            September, 1992
  ;       Major rewrite for improved speed   W. Landsman    April 2000
  ;       Accept up to 25 variables, use SCOPE_VARFETCH internally
  ;              W. Landsman   Feb 2010
  ;       Fix occasional integer overflow problem  V. Geers  Feb 2011
  ;       Quietly return if index is !null or undefined W.L. Aug 2011
  ;
  ;-
  On_error,2
  compile_opt idl2,strictarrsubs

  npar = N_params()
  nvar = npar-1
  if npar LT 2 then begin
    print,'Syntax - remove, index, v1, [v2, v3, v4,..., v25]'
    return
  endif

  if N_elements(index) EQ 0 then return

  vv = 'v' + strtrim( indgen(nvar)+1, 2)
  npts = N_elements(v1)

  max_index = max(index, MIN = min_index)

  if ( min_index LT 0 ) || (max_index GT npts-1) then message, $
    'ERROR - Index vector is out of range'

  if ( max_index Eq min_index ) then begin   ;Remove only 1 element?
    Ngood = 0
    if npts EQ 1 then message, $
      'ERROR - Cannot delete all elements from a vector'
  endif else begin


    ;  Begin case where more than 1 element is to be removed.   Use HISTOGRAM
    ;  to determine then indices to keep

    nhist = max_index - min_index +1

    hist = histogram( index)      ;Find unique index values to remove
    keep = where( hist EQ 0, Ngood ) + min_index

    if ngood EQ 0 then begin
      if ( npts LE nhist ) then message, $
        'ERROR - Cannot delete all elements from a vector'
    endif
  endelse

  imin = min_index - 1
  imax = max_index + 1
  i0 = (min_index EQ 0) + 2*(max_index EQ npts-1)
  case i0 of
    3: begin
      for i=0, nvar-1 do  $
        (SCOPE_VARFETCH(vv[i],LEVEL=0)) = $
        (SCOPE_VARFETCH(vv[i],LEVEL=0))[keep]
      return
    end

    1:  ii = Ngood EQ 0 ? imax + lindgen(npts-imax) : $
      [keep, imax + lindgen(npts-imax) ]
    2:  ii = Ngood EQ 0 ? lindgen(imin+1)               :  $
      [lindgen(imin+1), keep ]
    0:   ii = Ngood EQ 0 ? [lindgen(imin+1), imax + lindgen(npts-imax) ]  : $
      [lindgen(imin+1), keep, imax + lindgen(npts-imax) ]
  endcase

  for i=0,nvar-1 do  $
    (SCOPE_VARFETCH(vv[i],LEVEL=0)) =    $
    (SCOPE_VARFETCH(vv[i],LEVEL=0))[ii]

  return
end

Attachment.