IDL_原函数分析2

PRO read_ascat
  COMPILE_OPT idl2
  ;ENVI调用初始化
  ENVI,/restore_base_save_files
  ENVI_BATCH_INIT
  file=DIALOG_PICKFILE(path='D:\oilSpill\ASCAT\解压',TITLE='Select unc File', FILTER='*.nc')
  nid = NCDF_OPEN(file, /nowrite )
  id1 = NCDF_VARID(nid, 'wind_dir')
  NCDF_VARGET, nid, id1, wind_dir
  id2 = NCDF_VARID(nid, 'lat')
  NCDF_VARGET, nid, id2, lat
  id3 = NCDF_VARID(nid, 'lon')
  NCDF_VARGET, nid, id3, lon
  id4 = NCDF_VARID(nid, 'wind_speed')
  id5=NCDF_VARID(nid, 'time')
  NCDF_VARGET, nid, id5, time
  ;help,time
  NCDF_VARGET, nid, id4, wind_speed
  NCDF_ATTGET,nid,id1,'scale_factor',a1
  NCDF_ATTGET,nid,id2,'scale_factor',a2
  NCDF_ATTGET,nid,id3,'scale_factor',a3
  NCDF_ATTGET,nid,id4,'scale_factor',a4
  index1=WHERE(wind_speed EQ -32767)
  wind_speed[index1]=0.0
  index2=WHERE(wind_dir EQ -32767)
  wind_dir[index2]=0.0
  wind_dir =wind_dir*a1     ;风速向量x风速“影响因素”a1
  lat=lat*a2                ;纬度X纬度“影响因素”a2，下同理
  lon=lon*a3
   wind_speed=wind_speed*a4
  s=size(lon,/dimensions)   ;判断lon是几维数组，返回值  维度的个数。
  FOR i=0,s[0]-1 DO BEGIN   ;对一维度循环
    for j=0,s[1]-1 do begin ;二维度循环      加上上一条分析  就是循环这个二维数组所有值
      IF (lon[i,j] GT 180) then  begin  ;判断值是否大于180
     lon[i,j] = lon[i,j]-360          ;如果其中的值大于180  就 -360
    endif
  endfor
endfor
  index3=WHERE((lon GT 0.5) AND (lon LT 1))          ;得到大于0.5 小于1 的数组索引 ，下同理 ，利用这个”索引的值“，判断大于57.5的和小鱼58的值，再取到这个域的索引
  index4=WHERE((lat[index3] GT 57.5) AND (lat[index3] LT 58) )
  b=index3[index4]                        ;同理  利用上面的域  ，把域作为索引，取到索引对应的值
  print,b
  print,lon[b]
  print,lat[b]
  print,wind_speed[b]
  print,wind_dir[b]
  print,time[b]
  t=fltarr(2)
  sum1=TOTAL(wind_dir[index3[index4]])
  r1=where(wind_dir[index3[index4]] ne 0,count)
  c1=size(r1,/n_elements)
  dir_mean=FLOAT(sum1/c1)
  print,'研究区域像元数:',c1
  PRINT,'平均风向：',dir_mean
  sum2=TOTAL(wind_speed[index3[index4]])
  r2=WHERE(wind_speed[index3[index4]] NE 0,count)
  c2=SIZE(r2,/n_elements)
  spd_mean=FLOAT(sum2/c2)
  PRINT,'研究区域像元数:',c2
  PRINT,'平均风速：',spd_mean
  index5=WHERE((lon GT -1) AND (lon LT 3))
  index6=WHERE((lat[index5] GT 56) AND (lat[index5] LT 59) )
  b1=index5[index6]
  cgDisplay, 900, 700, Title='wind on Map'
  thisDevice = !D.Name
  cgMap_Set, /Cylindrical, /NoBorder, /NoErase, $
    Limit=[56,-1,59.5,3.5], Position=[0.1, 0.1, 0.75, 0.85]
  ;MAP_SET ,58,0,limit=[55,-10,65,10] ,CHARSIZE=0.8,/CYLINDRICAL ,E_GRID={LABEL:2}
  cgMap_Continents, Color=80B, /Fill
  spd_colors = ['purple', 'dodger blue', 'dark green',$
    'green', 'yellow']
  ;TVLCT, cgColor(spd_colors, /Triple), 1
  cgLoadCT,33
  TVLCT, cgColor('gray', /Triple), 0
  SET_PLOT,thisDevice
  ;symbol = cgSymcat(15)
  ;spd_colors = BYTSCL(wind_speed[b1],top=4)+1B
  spd_colors = BYTSCL(wind_speed[b1])
  ;cgplots,0.5,57.5, PSym=symbol,Color=255, SymSize=2.1
  r=size(b1,/n_elements)
  print,r
 ; FOR i=0,r-1 DO BEGIN
     ; IF ((wind_speed[b1[i]] NE 0 ))THEN BEGIN
       ; cgplots,lon[b1[i]], lat[b1[i]], PSym=symbol,Color=spd_colors[i], SymSize=1.4
     ; ENDIF
 ; ENDFOR
  ;FOR i=0,s[0]-1 DO BEGIN
 ;   FOR j=0,s[1]-1 DO BEGIN
 ; IF ((wind_speed[i,j] NE 0 )AND (lon[i,j] GT -10) AND (lon[i,j] LT 15) AND  (lat[i,j]  GT 50) AND (lat[i,j] LT 70)) THEN BEGIN
  ;   cgplots,lon[i,j], lat[i,j], PSym=symbol,Color=spd_colors[j*s[0]+i], SymSize=1.4
  ;ENDIF
;ENDFOR
;ENDFOR
  cgColorbar, /Vertical,Bottom=1, NColors=255, $
   Divisions=5, Minor=0, YTicklen=0.5, Range=[MIN(wind_speed[b1]),MAX(wind_speed[b1])], $
    /right, Title='m/s', Format='(F5.2)',Position=[0.84, 0.1, 0.87, 0.85]
    FOR i=0,r-1 DO BEGIN
      p=CONVERT_COORD([lon[b1[i]]],[lat[b1[i]]],/to_device)
      t[0]=ROUND(p[0])+25*SIN(wind_dir[b1[i]]*!pi/180)
      t[1]=ROUND(p[1])+25*COS(wind_dir[b1[i]]*!pi/180)
      IF ((wind_speed[b1[i]]) NE 0 )THEN BEGIN
        cgARROW, ROUND(p[0]),ROUND(p[1]), t[0], t[1],HTHICK=1.5, HSIZE=12.0,Color=spd_colors[i],/solid
      ENDIF
      i=i+1
  ENDFOR
  ;FOR i=0,s[0]-1 DO BEGIN
    ; for j=0,s[1]-1 do begin
    ; p=convert_coord( [lon[i,j]],[lat[i,j]],/to_device)
    ; t[0]=ROUND(p[0])+17*sin(wind_dir[i,j]*!pi/180)
     ;t[1]=ROUND(p[1])+17*cos(wind_dir[i,j]*!pi/180)
     ;if ((wind_speed[i,j] ne 0 )and (lon[i,j] GT -10) AND (lon[i,j] LT 15) and  (lat[i,j]  GT 50) AND (lat[i,j] LT 70)) then begin
    ;ARROW, ROUND(p[0]),ROUND(p[1]), t[0], t[1],HTHICK=0.05, HSIZE=7.0,/device
   ; endif
     ; j=j+2
    ; endfor
     ;i=i+1
    ;endfor
    file_name=file_basename(file)
    A=strmid(file_name,6,4)
    B=strmid(file_name,11,1)
    C=strmid(file_name,12,2)
    d=strmid(file_name,15,2)
    e=strmid(file_name,17,2)
    p1=CONVERT_COORD( [0.5],[57.5],/to_device)
    p2=CONVERT_COORD( [1.0],[58],/to_device)
    PLOTS,[p1[0],p1[0]],[p1[1],p2[1]],/DEVICE,color=255,thick=2
    PLOTS,[p1[0],p2[0]],[p1[1],p1[1]],/DEVICE,color=255,thick=2
    PLOTS,[p2[0],p2[0]],[p1[1],p2[1]],/DEVICE,color=255,thick=2
    PLOTS,[p1[0],p2[0]],[p2[1],p2[1]],/DEVICE,color=255,thick=2
   ;cgMap_Continents, Color='charcoal'
    cgMap_Grid, /Box, Color='charcoal',charsize=0.7
    cgText, 0.425, 0.925, /Normal, Alignment=0.5,'ASCAT '+A+'.'+b+'.'+c+' '+d+':'+e, Charsize=cgDefCharsize()*1.25
   copyimage=TVRD(true=1)
   write_jpeg,'D:\oilSpill\修改之后的结果图\ASCAT风场\ascat_new\'+b+c+d+e+'.JPG',copyimage,true=1
 end