| wiki | recipes | search |
| Title: | Gridding and Plotting Station Data |
| Author: | Luiz Rodrigo Tozzi <luizrodrigotozzi@gmail.com> |
| Created: | 2008-10-09 |
| Status: | WORK IN PROGRESS |
| Requires: |
Contents |
Problem
A common approach in Spatial Statiscs is to transform scattered and unsctructured data for the purpose of creating maps. Those maps can display data from various sources, like a weather station network or a georeferenced data mining procedure.
Meteorological station data are generated by different sources and in different formats and normally GrADS cannot access them automatically. Just like in previous recipes, the data needs to be converted to a proper format for ingestion in GrADS. We will use an objective analysis scheme to transform the station data to a regular grid so that we can contour it.
This recipe derives an operational procedure to grid the daily precipitation from a weather station network in Rio de Janeiro (Brazil), enabling its comparison to forecast fields on a single map.
Solution
The solution here was adapted from a question posted to the GRADSUSR list in January, 2008, by Qingfang Dai and it's solution (by JayaKrishnan PR). It uses a Fortran program to arrange the data and a C program to make it binary.
The following files will be used:
| Contains the daily precipitation data from all the stations in the network, for a given day. | |
| Contains any forecasted variable and, most important, the information of the grid we want the station data to be gridded to. | |
|
station2modelgrid.c |
Scripts and codes needed. |
Executing the main script
The shell file station2modelgrid.sh will be used to run all the steps automatically. The content of the script is bellow:
# station2modelgrid.sh # Luiz Rodrigo Tozzi - luizrodrigotozzi@gmail.com # stationdata.txt (unformatted) -> station2modelgrid.txt (formatted) gfortran -i ./station2modelgrid.f -o ./station2modelgridF ./station2modelgridF # station2modelgrid.txt (ASCII) -> station2modelgrid.bin (BINARY) gcc ./station2modelgrid.c -o ./station2modelgridC ./station2modelgridC # generate the product in GrADS stnmap -i ./station2modelgrid.ctl gradsc -blc "./station2modelgrid.gs"
Note that in the script I assume that all programs, including the Fortran and C compilers (gfortran and gcc, respectivelly), are in the path.
Formatting the data
The first thing the shell script will do is format the station data. The sample station data we provide is different from the one in 'Recipe-005: Creating GrADS binary station data from ASCII' because in the other recipe the data is annual and in this one the data can be hourly (with the restriction of one time per file). The file stationdata.txt looks like:
2007 12 01 00 station002 -22.992200 -43.232800 7
2007 12 01 00 station003 -22.955600 -43.166700 6
2007 12 01 00 station004 -22.985600 -43.245600 7
2007 12 01 00 station005 -22.931700 -43.221700 7
2007 12 01 00 station007 -22.985800 -43.188600 6
2007 12 01 00 station008 -22.923600 -43.267800 7
2007 12 01 00 station009 -22.817500 -43.210000 7
2007 12 01 00 station010 -22.843600 -43.274400 7
2007 12 01 00 station011 -22.872200 -43.338900 4
2007 12 01 00 station012 -22.821400 -43.336700 5
2007 12 01 00 station013 -22.886100 -43.465300 11
2007 12 01 00 station014 -22.891100 -43.308600 5
2007 12 01 00 station015 -22.910300 -43.366400 6
2007 12 01 00 station016 -22.897200 -43.194200 6
2007 12 01 00 station017 -22.972200 -43.223300 7
2007 12 01 00 station018 -22.998300 -43.302200 7
2007 12 01 00 station019 -22.945000 -43.362200 7
2007 12 01 00 station020 -22.975000 -43.411700 5
2007 12 01 00 station021 -23.052200 -43.577800 26
2007 12 01 00 station022 -22.866100 -43.583300 23
2007 12 01 00 station023 -22.912200 -43.684400 34
2007 12 01 00 station024 -22.889700 -43.278100 8
2007 12 01 00 station025 -22.826700 -43.403100 3
2007 12 01 00 station026 -23.014400 -43.518900 14
2007 12 01 00 station027 -22.903300 -43.561700 25
2007 12 01 00 station028 -22.968300 -43.711700 36
2007 12 01 00 station029 -22.950800 -43.235800 7
2007 12 01 00 station030 -22.823600 -43.523300 13
2007 12 01 00 station031 -23.009200 -43.424700 14
2007 12 01 00 station032 -22.937500 -43.187200 8
The code of station2modelgrid.f is bellow:
! PROGRAM: CreateStation2ModelGrid ! AUTHOR: Luiz Rodrigo Tozzi - luizrodrigotozzi@gmail.com ! DESCRIPTION: This fortran code converts a free-styled table with "station" values to a formatted table. This is the first step to create a GrADS station binary file. ! ! INPUT FILE LOOKS LIKE (separeted by any number of spaces): ! ! (...) ! yearwith4digits monthwith2digits daywith2digits hourwith2digits stationname latitude longitude value ! (...) ! program CreateStation2ModelGrid parameter (nl=9999) ! nl: MAXIMUM NUMBER OF STATIONS real lat(nl),lon(nl),var(nl) integer m2,y4,status character (LEN=3) m2name(12) character (LEN=2) d2,h2 character sid open(unit=1,file='stationdata.txt') open(unit=2,file='station2modelgrid.txt') open(unit=10,file='station2modelgrid.ctl') i=0 do read(1,*,END=21)y4,m2,d2,h2,sid,lat(i),lon(i),var(i) if((lat(i).eq.lon(i)).and.(lon(i).eq.0))exit i=i+1 20 end do 21 if(d2(1:1).eq." ") then d2(1:1)="0" endif if(d2(2:2).eq." ") then d2(2:2)=d2(1:1); d2(1:1)="0" endif if(h2(1:1).eq." ") then h2(1:1)="0" endif if(h2(2:2).eq." ") then h2(2:2)=h2(1:1); h2(1:1)="0" endif do 30 j=0,i ! Transforms the station name (sid) in an index number (i) if(m2.lt.10) then write(*,15)y4,"0",m2,d2,h2,j,lat(j),lon(j),var(j) write(2,15)y4,"0",m2,d2,h2,j,lat(j),lon(j),var(j) else write(*,14)y4,m2,d2,h2,j,lat(j),lon(j),var(j) write(2,14)y4,m2,d2,h2,j,lat(j),lon(j),var(j) endif 30 continue 14 format(i4,2x,i2,2x,a2,2x,a2,2x,i5,9x,f10.3,2x,f10.3,2x,f10.3) 15 format(i4,2x,a1,i1,2x,a2,2x,a2,2x,i5,9x,f10.3,2x,f10.3,2x,f10.3) !!!!! Generating CTL !! Picking the right m2 name m2name=(/'jan','feb','mar','apr','may','jun','jul', *'aug','sep','oct','nov','dec'/) write(10,*) "DSET ^station2modelgrid.bin" write(10,*) "DTYPE station " write(10,*) "STNMAP station2modelgrid.map" write(10,*) "ZDEF 1 1" write(10,*) "UNDEF -9.99e33" write(10,*) "TITLE Station Data" write(10,'(a16,a2,a1,a2,a3,i4,a5)') "TDEF 1 linear ",h2, *"z",d2,m2name(m2),y4," 12hr" write(10,*) "VARS 1" write(10,*) "var 0 99 **" write(10,*) "ENDVARS" stop end
This Fortran code creates the ASCII formatted file station2modelgrid.txt and a control file (station2modelgrid.ctl).
The only thing editable in the script, if necessary, is the maximum number of stations (lines of the file).
Transforming the data to a binary station file
The next step is to get the formatted ASCII data and transform it to binary. The station2modelgrid.c code is a really simple version of stngrads code in 'Recipe-005: Creating GrADS binary station data from ASCII'.
Here is the C code of station2modelgrid.c:
/* station2modelgrid.c Luiz Rodrigo Tozzi - luizrodrigotozzi@gmail.com */ #include <stdio.h> /* Structure that describes a report header in a stn file */ struct rpthdr { char id[8]; /* Station ID */ float lat; /* Latitude of Station */ float lon; /* Longitude of Station */ float t; /* Time in grid-relative units */ int nlev; /* Number of levels following */ int flag; /* Level independent var set flag */ } hdr; main () { FILE *ifile, *ofile; char rec[80]; int flag,year,month,day,hour,yrsav,mnsav,ddsav,hhsav,i; float val; /* Open files */ ifile = fopen ("station2modelgrid.txt","r"); ofile = fopen ("station2modelgrid.bin","wb"); if (ifile==NULL || ofile==NULL) { printf("Error opening files\n"); return; } /* Read, write loop */ flag = 1; while (fgets(rec,79,ifile)!=NULL) { /* Format conversion */ sscanf (rec,"%i %i %i %i",&year,&month,&day,&hour); sscanf (rec+26," %g %g %g",&hdr.lat,&hdr.lon,&val); for (i=0; i<8; i++) hdr.id[i] = rec[i+11]; /* Time group terminator if needed */ if (flag) { yrsav = year; mnsav = month; ddsav = day; hhsav = hour; flag = 0; } if (yrsav!=year || mnsav!=month || ddsav!=day) { hdr.nlev = 0; fwrite(&hdr,sizeof(struct rpthdr), 1, ofile); } yrsav = year; mnsav = month; ddsav = day; hhsav = hour; /* Write this report */ hdr.nlev = 1; hdr.flag = 1; hdr.t = 0.0; fwrite (&hdr,sizeof(struct rpthdr), 1, ofile); fwrite (&val,sizeof(float), 1, ofile); } hdr.nlev = 0; fwrite (&hdr,sizeof(struct rpthdr), 1, ofile); }
The code above creates station2modelgrid.bin, that can be read by GrADS through the station2modelgrid.ctl control file.
Gridding the station data
Once we created the binary station data, the station2modelgrid.gs GrADS script will be used to grid it using an objecive analysis scheme.
Here is the script:
* station2modelgrid.gs * Luiz Rodrigo Tozzi - luizrodrigotozzi@gmail.com **** Open the station file 'open station2modelgrid.ctl' 'query time' lin=sublin(result,1) time=subwrd(lin,3) 'set t 1' **** Open the numerical model file 'open wrf.ctl' 'set lon -43.75 -43.1' 'set lat -23.1 -22.7' 'set string 1' 'set xsize 750 550' 'set display color white' *1*'set mpdraw off' **** Check if the field has values > 0 'set gxout stat' 'd var' lin=sublin(result,8) say lin min=subwrd(lin,4) max=subwrd(lin,5) **** Plot the grid of the station value (precipitation) in the numerical model grid IF VALUES ARE > 0 'set gxout shaded' 'set grid off' 'set grads off' 'd oacres(slvl.2,var)' if(min=max) 'c' endif 'run cbarn' ****Plot the station value (precipitation) 'set digsiz 0.2' 'set line 1 1 6' *1*'shp_lines mun' 'set grads off' 'd var' **** Plot the mean sea level pressure from the numerical model grid 'set gxout contour' 'set grid off' 'set grads off' 'set cstyle 2' 'set ccolor 1' 'set cthick 1' 'set clab forced' 'd slvl.2' 'draw title Prec Stations+Pressure - 'time *2*'gxyat 'time'.png' 'printim 'time'.png' 'quit'
The script opens the station data file as first file and then opens the model data file. The objecive analysis occurs using the OACRES function and this script displays, in this order, the shaded plot of the gridded station data (daily precipitation total), the values of precipitation of each station in its real coordinate and a contour plot of mean sea level pressure.
Using the proper shapefile instead of GrADS maps and GxYat OpenGrADS User Defined eXtension instead of printim, the result of the script is this PNG file:
Discussion
The original motivation of this issue was to plot climatology station contours and grid it. Nowadays there are many other uses for this type of process, such as ploting data mining results against weather variables, in a georeferenced map.
The first part of this recipe (Formatting the data section) could be skipped if we could garantee the generation of the source station file (stationdata.txt, in this case) in the proper format. It's also possible to do this section using shell commands like grep and awk.
Remember that both control files (wrf.ctl and station2modelgrid.gs, in this case) must relate to the same time to be able to use the OACRES function.
See Also
OACRES GrADS function
Recipe-005: Creating GrADS binary station data from ASCII
No Warranty
Because the software provided in this Recipe is licensed free of charge, there is no warranty for it, to the extent permitted by applicable law. Except when otherwise stated in writing the authors and/or other parties provide the program "as is" without warranty of any kind, either expressed or implied, including, but not limited to, the implied warranties of merchantability and fitness for a particular purpose. The entire risk as to the quality and performance of the software is with you. Should the software prove defective, you assume the cost of all necessary servicing, repair or correction.
Copyright
This document has been placed in the public domain.
Last modified October 31, 2008 9:11 pm / Skin by Kevin Hughes
