|
Size: 6188
Comment:
|
← Revision 4 as of 2011-02-16 10:55:06 ⇥
Size: 1196
Comment:
|
| Deletions are marked like this. | Additions are marked like this. |
| Line 1: | Line 1: |
| == Code to read GSOD data for a specific station and convert it == |
|
| Line 4: | Line 2: |
| #!/usr/bin/env python #read GSOD data and store daily values for all stations in a single file #it is assumed that all files in the data directory belong to the same year! from pylab import * import glob import os |
#example how to make a sub-sampling of an array demsamp=10 #takes e.g. every etnth pixel ELEV=ELEV[::demsamp,::demsamp] }}} |
| Line 12: | Line 8: |
| inch2mm = 25.4 #conversion factor knots2ms = 0.51444 #convert from knots to m/s |
== Code to read GTOPO30 data == {{{#!python def read_gtopo30(filename): ''' routine to read GTOPO 30 data INPUT: filename without extension OUTPUT: DEM and LAT/LON GRID ''' print 'Reading GTOPO30 data: ', filename #/// read header /// hdrname = filename + '.HDR' data = [line.split() for line in open(hdrname).readlines()] hdr = dict(data) |
| Line 15: | Line 28: |
| year='2005' ddir='/net/nas2/export/eo/rawdata/station_data/GSOD/newdata/' + year + '/' odir='/net/nas2/export/eo/rawdata/station_data/GSOD/newdata/output/' |
ulx,uly = float(hdr['ULXMAP']),float(hdr['ULYMAP']) nx, ny = int(hdr['NCOLS']), int(hdr['NROWS']) dx, dy = float(hdr['XDIM']), float(hdr['YDIM']) nodata = float(hdr['NODATA']) |
| Line 19: | Line 33: |
| os.system('mkdir -p '+ odir) #!/// read GSOD station information /// f=open('all_stations.txt','r'); X=f.readlines() stat=[]; stat2=[]; name=[] ct1=[]; ct2=[]; ct3=[] lat=[]; lon=[]; elev=[] dummy=-99999. for i in range(len(X)): stat.append(X[i][0:6]) stat2.append(X[i][7:12]) name.append(X[i][13:42]) ct1=X[i][43:45]; ct2=X[i][46:48]; ct3=X[i][49:56] try: lat.append(float(X[i][57:64])/1000.) except: lat.append(dummy) try: lon.append(float(X[i][65:72])/1000.) except: lon.append(dummy) try: elev.append(float(X[i][73:79])) except: elev.append(dummy) lat=asarray(lat); lon=asarray(lon) lat[lat==dummy]=nan; lon[lon==dummy]=nan f.close() #/// process files /// #generate output arrays T=zeros((len(lat),366)) Tmax=zeros((len(lat),366)) Tmin=zeros((len(lat),366)) T[:,:]=nan; Tmax[:,:]=nan; Tmin[:,:]=nan YEAR=zeros((len(lat),366)); YEAR[:,:]=dummy MONTH=zeros((len(lat),366)) DAY=zeros((len(lat),366)) WDSP=zeros((len(lat),366)); WDSP[:,:]=dummy MXSPD=zeros((len(lat),366)); MXSPD[:,:]=dummy PRCP=zeros((len(lat),366)); PRCP[:,:]=dummy #/// process all files /// for k in range(len(lat)): #for k in range(100): s = stat[k] |
#/// read data /// elev = None file = open(filename + '.DEM','rb') a=array.array('h') a.read(file,nx*ny) a.byteswap() elev=reshape(a,(ny,nx)) elev=flipud(elev) a=None |
| Line 68: | Line 43: |
| print 'Processing station: ', k, '/', len(lat), s | #/// geometry /// dem_lon = arange(ulx,ulx+nx*dx,dx) dem_lat = arange(uly-ny*dy,uly,dy) LON,LAT=meshgrid(dem_lon,dem_lat) #elev[elev==nodata]= |
| Line 70: | Line 49: |
| istat = k files=glob.glob(ddir + s + '*-' + year + '.op') if len(files) != 1: if len(files) == 0: print 'Missing station ... ' + s continue else: print 'We have a problem here: ', files stop print 'Processing! ***' #/// process a single file /// file=files[0] f=open(file,'r') ls=f.readlines() for i in range(1,len(ls)): #for i in range(1,2): l=ls[i] sta =l[ 0:6 ] wban =l[ 7:12 ] year =l[14:18]; month =l[18:20]; day =l[20:22] t =l[26:30 ]; nt =l[31:33] #temperature [Fahrenheit] dew =l[37:41]; ndew =l[42:44] #dewpoint temperature (Fahrenheit) slp =l[46:52]; nslp =l[53:55] #sea level pressure stp =l[57:63]; nstp =l[64:66] #station pressure vis =l[67:73]; nvis =l[74:76] wdsp =l[78:83]; nwdsp =l[84:86] #wind speed in knots to tenth or meters per second mxspd =l[88:93] gust =l[95:100] tmax =l[102:108] #maximum temperature tmin =l[110:116] prcp =l[117:123] #precipitation [inch] - pflag =l[123:124] snow =l[125:130] if sta != s: print 'ERROR: station does not fit to data!', sta, s stop #/// julian day /// jd = datestr2num(year + '-' + month + '-' + day) - datestr2num(year + '-01-01') + 1 idx = int(jd-1) #/// store data /// if float(t) != 9999.9: T[istat,idx] = float(t) else: T[istat,idx] = nan if float(tmax) != 9999.9: Tmax[istat,idx] = float(tmax) else: Tmax[istat,idx] = nan if float(tmin) != 9999.9: Tmin[istat,idx] = float(tmin) else: Tmin[istat,idx] = nan YEAR[istat,idx] = int(year) MONTH[istat,idx] = int(month) DAY[istat,idx] = int(day) if float(wdsp) != 999.9: WDSP[istat,idx] = float(wdsp)*knots2ms #m/s else: WDSP[istat,idx] = nan if float(mxspd) != 999.9: MXSPD[istat,idx] = float(mxspd)*knots2ms else: MXSPD[istat,idx] = nan if float(prcp) != 99.99: PRCP[istat,idx] = float(prcp)*inch2mm #conversion to mm else: PRCP[istat,idx] = nan f.close() os.system('mv ' + file + ' ' + file + '.proc') #rename already processed files T=(T-32.) * 5. / 9. + 273.15 #temperature in degree Kelvin Tmax=(Tmax-32.) * 5. / 9. + 273.15#temperature in degree Kelvin Tmin=(Tmin-32.) * 5. / 9. + 273.15 #temperature in degree Kelvin T[isnan(T)] = dummy Tmax[isnan(Tmax)] = dummy Tmin[isnan(Tmin)] = dummy PRCP[isnan(PRCP)] = dummy WDSP[isnan(WDSP)] = dummy MXSPD[isnan(MXSPD)] = dummy #/// write output for each day /// for i in range(366): j=i+1 istr = str(j) if (j<100): istr = '0' + istr if (j<10): istr = '0' + istr print 'Writeing results for JD: ', istr f=open(odir + year + istr + '.out','w') #header 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') for k in range(len(lat)): if (~isnan(lon[k])) and (~isnan(lat[k]) and (YEAR[k,i] != dummy) ): 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' ) f.close() f.close() #now we can do the processing for each day ... #STN--- WBAN YEARMODA TEMP DEWP SLP STP VISIB WDSP MXSPD GUST MAX MIN PRCP SNDP FRSHTT #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 |
return LON,LAT,elev |
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