Code to read GSOD data for a specific station and convert it

   1 #!/usr/bin/env python
   2 #read GSOD data and store daily values for all stations in a single file
   3 #it is assumed that all files in the data directory belong to the same year!
   4 from pylab import *
   5 import glob
   6 import os
   7 
   8 
   9 inch2mm = 25.4 #conversion factor
  10 knots2ms = 0.51444 #convert from knots to m/s
  11 
  12 year='2005'
  13 ddir='/net/nas2/export/eo/rawdata/station_data/GSOD/newdata/' + year + '/'
  14 odir='/net/nas2/export/eo/rawdata/station_data/GSOD/newdata/output/'
  15 
  16 os.system('mkdir -p '+ odir)
  17 
  18 #!/// read GSOD station information ///
  19 f=open('all_stations.txt','r'); X=f.readlines()
  20 stat=[]; stat2=[]; name=[]
  21 ct1=[]; ct2=[]; ct3=[]
  22 lat=[]; lon=[]; elev=[]
  23 dummy=-99999.
  24 for i in range(len(X)):
  25     stat.append(X[i][0:6])
  26     stat2.append(X[i][7:12])
  27     name.append(X[i][13:42])
  28     ct1=X[i][43:45]; ct2=X[i][46:48]; ct3=X[i][49:56]
  29     try:
  30         lat.append(float(X[i][57:64])/1000.)
  31     except:
  32         lat.append(dummy)
  33     try:
  34         lon.append(float(X[i][65:72])/1000.)
  35     except:
  36         lon.append(dummy)
  37     try:
  38         elev.append(float(X[i][73:79]))
  39     except:
  40         elev.append(dummy)
  41 lat=asarray(lat); lon=asarray(lon)
  42 lat[lat==dummy]=nan; lon[lon==dummy]=nan
  43 f.close()
  44 
  45 #/// process files ///
  46 
  47 #generate output arrays
  48 T=zeros((len(lat),366))
  49 Tmax=zeros((len(lat),366))
  50 Tmin=zeros((len(lat),366))
  51 T[:,:]=nan; Tmax[:,:]=nan; Tmin[:,:]=nan
  52 YEAR=zeros((len(lat),366)); YEAR[:,:]=dummy
  53 MONTH=zeros((len(lat),366))
  54 DAY=zeros((len(lat),366))
  55 WDSP=zeros((len(lat),366)); WDSP[:,:]=dummy
  56 MXSPD=zeros((len(lat),366)); MXSPD[:,:]=dummy
  57 PRCP=zeros((len(lat),366)); PRCP[:,:]=dummy
  58 
  59 
  60 #/// process all files ///
  61 for k in range(len(lat)):
  62 #for k in range(100):
  63     s = stat[k]
  64     
  65     print 'Processing station: ', k, '/', len(lat), s
  66     
  67     istat = k
  68     files=glob.glob(ddir + s + '*-' + year + '.op')
  69     if len(files) != 1:
  70         if len(files) == 0:
  71             print 'Missing station ... ' + s
  72             continue
  73         else:
  74             print 'We have a problem here: ', files
  75             stop
  76     print 'Processing! ***'
  77 
  78     #/// process a single file ///
  79     file=files[0]
  80     f=open(file,'r')
  81     ls=f.readlines()
  82     for i in range(1,len(ls)):
  83     #for i in range(1,2):
  84         l=ls[i]
  85         sta   =l[ 0:6 ]
  86         wban  =l[ 7:12 ]
  87         year  =l[14:18];  month =l[18:20]; day   =l[20:22]
  88         t     =l[26:30 ]; nt    =l[31:33] #temperature [Fahrenheit]
  89         dew   =l[37:41];  ndew  =l[42:44] #dewpoint temperature (Fahrenheit)
  90         slp   =l[46:52];  nslp  =l[53:55] #sea level pressure
  91         stp   =l[57:63];  nstp  =l[64:66] #station pressure
  92         vis   =l[67:73];  nvis  =l[74:76]
  93         wdsp  =l[78:83];  nwdsp =l[84:86] #wind speed in knots to tenth or meters per second
  94         mxspd =l[88:93]
  95         gust  =l[95:100]
  96         tmax  =l[102:108] #maximum temperature
  97         tmin  =l[110:116]
  98         prcp  =l[117:123] #precipitation [inch] -
  99         pflag =l[123:124]
 100         snow  =l[125:130]
 101         
 102         if sta != s:
 103             print 'ERROR: station does not fit to data!', sta, s
 104             stop
 105         
 106         #/// julian day ///
 107         jd = datestr2num(year + '-' + month + '-' + day) - datestr2num(year + '-01-01') + 1
 108         
 109         idx = int(jd-1)
 110         
 111         #/// store data ///
 112         if float(t) != 9999.9:
 113             T[istat,idx]      = float(t) 
 114         else:
 115             T[istat,idx] = nan
 116         if float(tmax) != 9999.9:
 117             Tmax[istat,idx]   = float(tmax) 
 118         else:
 119             Tmax[istat,idx]   = nan
 120         if float(tmin) != 9999.9:
 121             Tmin[istat,idx]   = float(tmin)
 122         else:
 123             Tmin[istat,idx]   = nan
 124             
 125         YEAR[istat,idx]   = int(year)
 126         MONTH[istat,idx]  = int(month)
 127         DAY[istat,idx]    = int(day)
 128         
 129         if float(wdsp) != 999.9:
 130             WDSP[istat,idx]   = float(wdsp)*knots2ms #m/s
 131         else:
 132             WDSP[istat,idx] = nan
 133             
 134         if float(mxspd) != 999.9:
 135             MXSPD[istat,idx]  = float(mxspd)*knots2ms
 136         else:
 137             MXSPD[istat,idx] = nan
 138             
 139         if float(prcp) != 99.99:
 140             PRCP[istat,idx]   = float(prcp)*inch2mm #conversion to mm
 141         else:
 142             PRCP[istat,idx] = nan
 143         
 144     f.close()
 145     os.system('mv ' + file + ' ' + file + '.proc') #rename already processed files
 146     
 147     
 148 T=(T-32.) * 5. / 9. + 273.15      #temperature in degree Kelvin 
 149 Tmax=(Tmax-32.) * 5. / 9. + 273.15#temperature in degree Kelvin
 150 Tmin=(Tmin-32.) * 5. / 9. + 273.15 #temperature in degree Kelvin
 151 
 152 T[isnan(T)]       = dummy
 153 Tmax[isnan(Tmax)] = dummy
 154 Tmin[isnan(Tmin)] = dummy
 155 PRCP[isnan(PRCP)] = dummy
 156 WDSP[isnan(WDSP)] = dummy
 157 MXSPD[isnan(MXSPD)] = dummy
 158 
 159 
 160 #/// write output for each day ///
 161 for i in range(366):
 162     j=i+1
 163     istr = str(j)
 164     if (j<100):
 165         istr = '0' + istr
 166     if (j<10):
 167         istr = '0' + istr
 168         
 169     print 'Writeing results for JD: ', istr
 170     f=open(odir + year + istr + '.out','w') 
 171     #header
 172     f.write('YEAR \t MONTH \t DAY \t LON \t LAT \t ELEV \t Tmean_K \t Tmin_K \t Tmax_K \t WINDMEAN_MS \t WINDMAX_MS \t PRECIP_MM \n')
 173     for k in range(len(lat)):
 174         if (~isnan(lon[k])) and (~isnan(lat[k]) and (YEAR[k,i] != dummy)     ):
 175             f.write(str(int(YEAR[k,i])) + '\t' + str(int(MONTH[k,i])) + '\t' + str(int(DAY[k,i])) + '\t' + str(lon[k]) + '\t' +  str(lat[k]) + '\t' + str(elev[k]) + '\t' +  str(T[k,i]) + '\t' +  str(Tmin[k,i]) + '\t' +  str(Tmax[k,i]) + '\t' + str(WDSP[k,i]) + '\t' + str(MXSPD[k,i]) + '\t' +  str(PRCP[k,i]) + '\n' )
 176     f.close()
 177 
 178 f.close()
 179 
 180 
 181 #now we can do the processing for each day ...
 182 
 183 #STN--- WBAN   YEARMODA    TEMP       DEWP      SLP        STP       VISIB      WDSP     MXSPD   GUST    MAX     MIN   PRCP   SNDP   FRSHTT
 184 #276070 99999  20050101    18.3 24    14.5 24  9999.9  0  9999.9  0    4.6 24   12.9 24   17.5   25.3    26.6*    8.6* 99.99  999.9  011100

Code to read GTOPO30 data

   1 def read_gtopo30(filename):
   2     '''
   3     routine to read GTOPO 30 data
   4     
   5     INPUT:
   6     filename without extension
   7     
   8     OUTPUT:
   9     DEM and LAT/LON GRID
  10     '''
  11     
  12     print 'Reading GTOPO30 data: ', filename
  13     
  14     #/// read header ///
  15     hdrname = filename + '.HDR'
  16     data = [line.split() for line in open(hdrname).readlines()]
  17     hdr = dict(data)
  18 
  19     ulx,uly = float(hdr['ULXMAP']),float(hdr['ULYMAP'])
  20     nx, ny = int(hdr['NCOLS']), int(hdr['NROWS'])
  21     dx, dy = float(hdr['XDIM']), float(hdr['YDIM'])
  22     nodata = float(hdr['NODATA'])
  23 
  24     #/// read data ///
  25     elev = None
  26     file = open(filename + '.DEM','rb')
  27     a=array.array('h')
  28     a.read(file,nx*ny)
  29     a.byteswap()
  30     elev=reshape(a,(ny,nx))
  31     elev=flipud(elev)
  32     a=None
  33     
  34     #/// geometry ///
  35     dem_lon = arange(ulx,ulx+nx*dx,dx)
  36     dem_lat = arange(uly-ny*dy,uly,dy)
  37     LON,LAT=meshgrid(dem_lon,dem_lat)
  38     #elev[elev==nodata]=
  39     
  40     return LON,LAT,elev