Friday, November 20, 2015

ENVI+IDL: 批量提取多幅影像共同区域的方法

问题：已知多幅裁剪之后的影像，他们存在着多一行（或列）的问题，将其中行（列）最小的影像作为多幅影像的共同区域，达到批量提取其他影像的目的？
方法：ENVI+IDL。
练习数据中，mask.tif是行（列）最小的影像，即为多幅影像的共同区域。1比mask.tif多一行（列），代码运行的结果是result，行（列）与mask.tif一致，其第一波段即是mask.tif，其他波段为1影像。
文件test.pro是代码的入口。代码运行结果与ENVI操作结果一致。
注意不同文件之间数据类型差别，可能在批处理过程中要适当修改代码中选择数据类型的规则（out_dt=max(in_dt)）。Referencing: IDL Data Types.
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
;;Created by LI Xu
;;Version 1.0
;;18 October, 2014

pro CommonRange_Batch

  print, 'Time Begin:'+' '+systime()
  begintime=SYSTIME(1)

  ;;Source Directory
  SouDir='E:\Tools\rad\Rad'
  ;;Destination Directory
  DesDir='E:\Tools\rad\NEW'
  ;;One Image
  imgpath='E:\Tools\rad\pro\mask.tif'
  
  
  ;;Retrieve all files
  filespath=file_search(SouDir, '*.tif', count=num_files, /test_regular)
  for ii=0, num_files-1 do begin
    file=filespath(ii)
    ;;file=strsplit(file, '.', /extract)
    ;;file=file(0)
    
    filename=strsplit(file, '\', /extract, count=strcount)
    filename=filename(strcount-1)
    otImage=DesDir+'\'+filename
    inImage=strarr(2)
    inImage[0]=imgpath
    inImage[1]=file
    
    CommonRange, inImage, otImage
    
    
  endfor
  

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




end


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
;+
; :Author: LiXu
;-

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

;;Description:
;;Make files in same Column and Row,
;;and upper left coodinate.
;;Reference:
;;http://www.exelisvis.com/Support/HelpArticlesDetail/TabId/219/ArtMID/900/ArticleID/4572/4572.aspx

pro CommonRange, imgpath, otimg

  ;;Debug Bolck
;  imgpath=strarr(2)
;  imgpath[0]='C:\AAAAA\practice\mask.tif'
;  imgpath[1]='C:\AAAAA\practice\1'
;  otimg='C:\AAAAA\practice\result'
  
  
  
  compile_opt idl2
  envi, /restore_base_save_files
  
  ;;Number of files
  num_files=N_Elements(imgpath)
  ;;open and gather file information in arrays
  in_fid=lonarr(num_files)
  in_nb=lonarr(num_files)
  in_dims=lonarr(5, num_files)
  in_dt=lonarr(num_files)
  
  for i=0L, num_files-1 do begin
    envi_open_file, imgpath[i], r_fid=r_fid
    if (r_fid eq -1) then begin
      print, imgpath[i]+' Open failed!'
      return
    endif

    envi_file_query, r_fid, ns=ns, nl=nl, nb=nb, dims=dims, data_type=dt
    in_fid[i]=r_fid
    in_nb[i]=nb
    in_dims[*,i]=dims
    in_dt[i]=dt
  endfor
  
  
  ;set up output fid, pos, and dims arrays
  out_fid = replicate(in_fid[0], in_nb[0])
  for i=1, num_files-1 do out_fid = [out_fid, replicate(in_fid[i], in_nb[i])]

  out_pos = lindgen(in_nb[0])
  for i = 1, num_files-1 do out_pos = [out_pos, lindgen(in_nb[i])]

  rep_dims = (intarr(in_nb[0])+1) # in_dims[*, 0]
  for i = 1, num_files-1 do $
    rep_dims = [rep_dims, (intarr(in_nb[i]) + 1) # in_dims[*, i]]
  out_dims = transpose(rep_dims)
  
  ;set the output projection and pixel size from the first file.
  ;save the result to disk and use max data type
  out_proj = envi_get_projection(fid=in_fid[0], pixel_size=out_ps)
  out_dt = min(in_dt)
  out_name=otimg
  
  ;call the layer stacking routine.
  ;Use nearest neighbor for the interpolation method.
  envi_doit, 'envi_layer_stacking_doit', fid=out_fid, pos=out_pos, dims=out_dims, $
    out_dt=out_dt, out_name=out_name, interp=0, out_ps=out_ps, $
    out_proj=out_proj, r_fid=r_fid,  /EXCLUSIVE

 

  
  
  
  
end
Attachment.
Posted by at No comments:
Labels: ,

Thursday, November 19, 2015

Python+GDAL: 文件格式的转换

空间数据的格式千变万化，常见的有ENVI Standard、GEOTIFF、Erads Image等等（栅格文件列表，矢量文件列表），以下应用GDAL工具转换文件的格式，省时高效。
举个栗子，ENVI Standard转换为GEOTIFF，参考How to call gdal_translate from Python code? Answered by Max。

Note

os.system (command): Execute the command (a string) in a subshell. This is implemented by calling the Standard C function system(). On Windows, the return value is that returned by the system shell after running command. 实际上，示例是Python调用C+GDAL的代码。
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
##Created by LI Xu
##Version 1.0
##7 October, 2014


##Convert ENVI Standard to GEOTIFF

##http://gis.stackexchange.com/questions/42584/how-to-call-gdal-translate-from-python-codearray-returns-just-nan-values-when-trying-to-read-envi-file


import gdal
import ogr
import os
import datetime
import time





def IsSubString(SubStrList,Str):  
 
    flag=True  
    for substr in SubStrList:  
        if not(substr in Str):  
            flag=False  
  
    return flag  
#~ #----------------------------------------------------------------------


def GetFileList(FindPath,FlagStr=[]):  

    import os  
    FileList=[]  
    FileNames=os.listdir(FindPath)  
    if (len(FileNames)>0):  
       for fn in FileNames:  
           if (len(FlagStr)>0):  
               
               if (IsSubString(FlagStr,fn)):  
                   fullfilename=os.path.join(FindPath,fn)  
                   FileList.append(fullfilename)  
           else:  
               
               fullfilename=os.path.join(FindPath,fn)  
               FileList.append(fullfilename)  
  
    
    if (len(FileList)>0):  
        FileList.sort()  
  
    return FileList  

def convert_format(infile, otfile, oFormat):
    os.system("gdal_translate -of" +" " +oFormat+" "+ infile + " " +  otfile)

##Main
begintime=time.strftime("%Y-%m-%d %H:%M:%S")
print 'Time Begin:'+begintime
starttime = datetime.datetime.now()

#Source Directory
SouDir=r'E:\Tools\rad\NEW'
#Destination Directory
DesDir=r'E:\Tools\rad\Rad'


#Retrieve files
FlagStr='.hdr'
files=GetFileList(SouDir,FlagStr)
for file in files:
    #print file
    str_file=file.split('.')
    infile=str_file[0]+'.tif'
    #print infile
    oformat="GTiff"
    ##Output Image
    filename=infile.split('\\')
    filename=filename[len(filename)-1]
    filename=filename.split('.')
    otimg=DesDir+'\\'+filename[0]+'.tif'
    print otimg
    convert_format(infile, otimg, oformat)
    



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'
Attachment.
Matlab+GDAL版本.
Posted by at No comments:
Labels: , ,

ENVI+IDL: Resampling

问题：已知一个影像，对另一影像像元分辨率按照已知影像进行重采样？
方法：ENVI+IDL。
练习数据包括两张GEOTIFF文件，具有相同的空间参考，STANDARD.tif是已知影像，像元分辨率：1000*1000METERS；reprojected.tif是被重采样影像，像元分辨率：799.232051*799.232051METERS。
（原来，我尝试Python+GDAL进行重投影和重采样，但没发现Python+GDAL重采样的代码，所以就分别做处理了。）
reprojected.tif与resampled.tif前后转换对比，如图 1（似乎没什么不同？）。
图 1
注意：土地利用数据的重采样方法应选择最邻近法，它基于离散数据在插值过程中不会产生新的数值，适用于类似土地利用分类数据的插值操作（NEAREST—Performs a nearest neighbor assignment, is the fastest of the interpolation methods. It is used primarily for discrete data, such as a land-use classification, since it will not change the values of the cells. The maximum spatial error will be one-half the cell size.）。 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
;;Created by LI Xu
;;Version 1.0
;;16 october, 2014


pro Resample_Doit_Batch
  print, 'Time Begin:'+' '+systime()
  begintime=SYSTIME(1)

  ;;Standard Image
  stdimgpath='E:\Tools\rad\pro\mask.tif'
  ;;Source Directory
  SouDir='E:\Tools\rad\Rad'
  ;;Destination Directory
  DesDir='E:\Tools\rad\NEW'
  ;;Interpolation method
  ;;Interp
  ;;0: Nearest neighbor
  ;;1: Bilinear
  ;;2: Cubic convolution
  ;;3: Pixel aggregate
  intermed=3
  
  ;;Retrieve all files
  files=file_search(SouDir, '*.tif', count=file_counts, /test_regular)
  for ii=0, file_counts-1 do begin
    file=files(ii)
    filename=strsplit(file, '\', count=str_num ,/EXTRACT)
    filename=filename(str_num-1)
    otfile=DesDir+'\'+filename
    ;print, otfile
    Resample_Doit, stdimgpath, file, intermed, otfile
  endfor
  
  print, 'End!'
  print, 'Time End:'+' '+systime()
  endtime=systime(1)
  timespan=endtime-begintime
  print, 'Time Span:'+' '+string(timespan)+' s'

end


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
;;Created by LI Xu
;;Version 1.0
;;16 October, 2014

;;Resampling x/y resolution size
;;Reference:
;;http://www.cnblogs.com/myyouthlife/archive/2012/07/01/2571785.html

;;Modified by LI Xu
;;Version 1.1
;;October 19, 2015
;;Add the variable for interpolating method


pro Resample_Doit, standardtif, resampltif, intemed, ottif
  
  ;;standardtif='D:\ForLIU\MOD_NDVI.tif'
  ;;resampltif='D:\ForLIU\LI_reprojected.tif'
  ;;ottif='D:\ForLIU\LI_resampled.tif'
  
  ;;*********************************
  compile_opt IDL2
  envi, /restore_base_save_files
  
  ;;Open the tifs
  ;;MODELPIXELSCALETAG DOUBLE[3] XYZ
  stdimg=read_tiff(standardtif, GEOTIFF=GEOVAR)
  tarXsize=GEOVAR.MODELPIXELSCALETAG(0)
  tarYsize=GEOVAR.MODELPIXELSCALETAG(1)
  envi_open_file, resampltif, r_fid=fid_resampled
  resampltif=read_tiff(resampltif, GEOTIFF=GEOVAR)
  nowXsize=GEOVAR.MODELPIXELSCALETAG(0)
  nowYsize=GEOVAR.MODELPIXELSCALETAG(1)
  envi_file_query, fid_resampled, dims=dims, nb=nb
  ;;DIMS[0]: A pointer to an open ROI; use only in cases where ROIs define the spatial subset. Otherwise, set to -1L.
  ;;DIMS[1]: The starting sample number. The first x pixel is 0.
  ;;DIMS[2]: The ending sample number
  ;;DIMS[3]: The starting line number. The first y pixel is 0.
  ;;DIMS[4]: The ending line number
  pos = lindgen(nb)
  
  Xsize=tarXsize/nowXsize
  Ysize=tarYsize/nowYsize
  envi_doit, 'resize_doit', $
    fid=fid_resampled, pos=pos, dims=dims, $
    interp=intemed, rfact=[Xsize, Ysize], $
    out_name=ottif, r_fid=r_fid
  ;;Interp  
  ;;0: Nearest neighbor
  ;;1: Bilinear
  ;;2: Cubic convolution
  ;;3: Pixel aggregate

end

References

[1] How does the Pixel Aggregate method work when resizing data with ENVI ?
Posted by at No comments:
Labels: ,

Tool: Tutorials Point

Summary

Tutorials Point originated from the idea that there exists a class of readers who respond better to online content and prefer to learn new skills at their own pace from the comforts of their drawing rooms.
Fig. 1
在Fig. 1当中的搜索框框中键入语言类型，如Matlab，随后选择Execute Matlab/Octave Online，跳跃至新页面就可以输入代码并在线测试代码，见Fig. 2。
Fig. 2
除Matlab之外，该网站还包括Python、C#、Java等语言编译器，在线编译器的缺点也很明显，速度较慢，因而比较适合临时编译代码。

References

Posted by at No comments:
Labels: ,

Python+GDAL: Reprojection

问题：已知一个影像，将另一影像按已知影像进行重投影？
方法：Python+GDAL。
练习数据包括两张GEOTIFF文件，STANDARD.tif是已知影像，空间参考信息（Spatial Reference）：PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]；tested.tif是另一影像，转换前的空间参考信息：GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]。
reprojected.tif是转换结果，空间参考信息：PROJCS["WGS 84 /UTM zone 50N",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]。
tested.tif与reprojected.tif前后转换对比，如图 1。
图 1
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
##Created by LI Xu
##Version 1.0
##16 October, 2014

##Discrption:
##Reproject a projected Image to other from a single Image
##Reference:
##http://pcjericks.github.io/py-gdalogr-cookbook/projection.html
##https://svn.osgeo.org/gdal/trunk/autotest/alg/reproject.py

import gdal
import ogr
import osr
import os
import datetime
import time
begintime=time.strftime("%Y-%m-%d %H:%M:%S")
print 'Time Begin:'+begintime
starttime = datetime.datetime.now()


gdal.AllRegister()

def reproject_img(stdimg, reprojimg, otimg):
    ##Retrieve the Spatial Reference of Standard Image
    standard=gdal.Open(stdimg)
    pszTarWKT=standard.GetProjection()
    print pszTarWKT

    ##Retrieve the Spatial Reference of Target Image
    tested=gdal.Open(reprojimg)
    pszSouWKT=tested.GetProjection()
    print pszSouWKT

    reproj_file = gdal.AutoCreateWarpedVRT(tested, pszSouWKT, pszTarWKT)
    reprojected=gdal.ReprojectImage(tested, reproj_file, pszSouWKT, pszTarWKT)
    reproj_attributes = reproj_file.GetGeoTransform()

    driver = gdal.GetDriverByName("GTiff")
    dest_file = driver.CreateCopy(otimg, reproj_file, 0)
    print dest_file.GetProjection()

    ##Close
    standard=None
    tested=None
    dest_file=None

    print otimg+' done!'

##Main
#Standard Image
StdImage=r'D:\NEW\STANDARD.tif'
#Source Directory
SouDir=r'D:\SPOT_Mon_Year'
#Destination Directory
DesDir=r'D:\NEW'

#Retrieve all tiffs
files=os.listdir(SouDir)
for file in files:
    filepath=SouDir+'\\'+file
    otfilepath=DesDir+'\\'+'SPOT_'+file
    reproject_img(StdImage, filepath,otfilepath)
    

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'
Attachment.
Matlab+GDAL版本.The 2nd Version 20230924.The 3rd Version 20240715.
Python批处理代码.
How to set up Python+GDAL? 2014-06-16
Posted by at No comments:
Labels: ,

Wednesday, November 18, 2015

Geography: NECT

Summary

The Northeast China Transect (NECT), one of the mid-latitude IGBP terrestrial transects, runs in parallel to 43º30'N and ranges from 42º to 46ºN and from 106º to 134ºE. The major global change gradient is precipitation ranging from 600-1000 mm in the east, 300-600 mm in the middle, and 100-300 mm in the west. Due to the steep moisture gradient, vegetation along the transect varies gradually from temperature evergreen conifer-deciduous broad leaf mixed forests, decidous broad leaf forests, woodlands, and shrublands in the east to typical steppes and desert steppes in the west, with agricultural fields, temperature savannas and meadow steppes in the middle. A secondary driving gradient is land use intensity from forest regions in the east, to agriculture in the middle, to pastoral areas in the west. The transect is an important base of forestry, agriculture and pastoral productions in China, producing wood, hay, grain crops (maize, soybean, wheat and rice) and cattle (leather, wool and milk).
Attachment.

References

[1] Ni, J., & Zhang, X.-S. (2000). Climate variability, ecological gradient and the Northeast China Transect (NECT). Journal of Arid Environments, 46(3), 313-325. doi: http://dx.doi.org/10.1006/jare.2000.0667.
[2] M_Map: A mapping package for Matlab.
[3] ArcGIS: Create polygon using file containing point coordinates.
Posted by at No comments:
Labels:

Monday, November 16, 2015

Matlab: Conversion from Arc/Info Grid to GeoTiff Format Program

Summary

记录ANUSPLIN插值输出结果在投影坐标系下的ASCII文件转换为GeoTiff格式文件。在前面的帖文中，我已经介绍Arc/Info Grid（Geographic Coordinate System）转换为Geotiff代码，此处追加Arc/Info Grid（Projected Coordinate System）转换为Geotiff代码。

Note

代码要求必须有一个Geotiff文件，它包含待转换Arc/Info Grid文件的全部地理信息，这些信息将在代码中赋给转换的矩阵进而输出为Geotiff文件。
Fig. 1
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
% Created by LI Xu
% Version 1.0
% October 26, 2015

% Description
% Main program for conversion to geotiff format 
% from Arc/Info Grid

% 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';
[~, geo]=geotiffread(maskpath);
info=geotiffinfo(maskpath);

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

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

for ii=1:length(files)
    filename=files(ii).name;
    % If exist
    otname=strsplit(filename, '.');
    otname=otname{1};
    otpath=[DesDir, '/', otname, '.tif'];
    
    if exist(otpath, 'file')
        disp([num2str(ii), ': ', otname, '.tif']);
        continue;
    end
    
    
    filepath=[SouDir, '/', filename];
    values=GetValues(filepath);
    geoattrs=GetProperty(filepath);
    inmat=reshape(values, geoattrs.NCOLS, geoattrs.NROWS);
    inmat=inmat';
    
    % Export
    
    geotiffwrite(otpath, inmat, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp([num2str(ii), ': ', otname, '.tif']);
end

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


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
% 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


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
% 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
Attachment.
Posted by at No comments:
Labels: ,
Subscribe to: Posts (Atom)