/* * $Id: st04c.c,v 1.5 1997/04/09 21:05:09 haley Exp $ */ /*********************************************************************** * * * Copyright (C) 1996 * * University Corporation for Atmospheric Research * * All Rights Reserved * * * ***********************************************************************/ /* * File: st04c.c * * Author: David Brown (converted by Mary Haley) * National Center for Atmospheric Research * PO 3000, Boulder, Colorado * * Date: Thu June 27 9:47:28 MST 1996 * * Description: This example shows a StreamlinePlot of 500 mb wind * vector data overlaid on a MapPlot. The streamlines * are drawn over a VectorPlot of surface winds colored * by surface pressure that in turn is drawn over a filled * ContourPlot of surface temperature. Different intervals * of the "temp1" colormap are used to color the contour * levels and the vectors. * The data represents 15 days of weather over North * America in January, 1996. * The data is extracted from NMC forcast data produced * at 12 hour intervals and converted to netcdf format * by Unidata. Most of the time steps in the files * extracted from the original data are taken from the * 0 and 6 hour forecast times. However, because some of * the original files were lost, certain time steps come * from longer range forcasts. Also, several steps had * to be excluded from the frame set because the data is * defective. The result is that there is an * apparent discontinuity between some of the frames * when the output is animated. */ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /* * Depending on the value of the TIMESTEPS variable declared below, * this example example can generate up to 61 frames from the 64 * timesteps in the data files. As shipped, only the first 20 frames * are created. To see the complete plot uncomment the second * assignment to TIMESTEPS. Some systems may not have enough physical * memory to allow all frames to be viewed as an animation. * * #define TIMESTEPS 64 */ #define TIMESTEPS 20 /* * Initialize netCDF filenames */ char *cdffiles[6] = {"Ustorm.cdf","Vstorm.cdf","Pstorm.cdf","Tstorm.cdf","U500storm.cdf","V500storm.cdf"}; main(int argc, char *argv[]) { int NCGM=1, X11=0, PS=0; int i, j, time, d, h; int appid, wid, cnid, vcid, stid, txid, amid, mpid, tmid, stdmid; long stid_len; int vfield, vfield2, sfield, sfield2; int rlist, len_dims[2]; long strt[1], cnt[1]; long latlen, lonlen; long timelen; int *timestep; int ncid[6], uid, vid, u5id, v5id, pid, tid; int latid, lonid; float *lon, *lat; float *X, *Y; char filename[256]; char *rftime; const char *dir = _NGGetNCARGEnv("data"); char hour[3], day[3], mainstring[17]; extern void get_2d_array(float *, long, long, int, int, long); /* * Initialize the high level utility library */ NhlInitialize(); /* * Create an application object. */ rlist = NhlRLCreate(NhlSETRL); NhlRLClear(rlist); NhlRLSetString(rlist,NhlNappUsrDir,"./"); NhlRLSetString(rlist,NhlNappDefaultParent,"True"); NhlCreate(&appid,"st04",NhlappClass,NhlDEFAULT_APP,rlist); if (NCGM == 1) { /* * Create a meta file workstation. */ NhlRLClear(rlist); NhlRLSetString(rlist,NhlNwkMetaName,"./st04c.ncgm"); NhlRLSetString(rlist,NhlNwkColorMap,"temp1"); NhlCreate(&wid,"st04Work", NhlncgmWorkstationClass,NhlDEFAULT_APP,rlist); } else if (X11) { /* * Create an X workstation. */ NhlRLClear(rlist); NhlRLSetInteger(rlist,NhlNwkPause,True); NhlRLSetString(rlist,NhlNwkColorMap,"temp1"); NhlCreate(&wid,"st04Work",NhlxWorkstationClass,appid,rlist); } else if (PS) { /* * Create a PS workstation. */ NhlRLClear(rlist); NhlRLSetString(rlist,NhlNwkPSFileName,"st04c.ps"); NhlRLSetString(rlist,NhlNwkColorMap,"temp1"); NhlCreate(&wid,"st04Work",NhlpsWorkstationClass,appid,rlist); } /* * Open the netCDF files. */ for( i = 0; i <= 5; i++ ) { sprintf( filename, "%s/cdf/%s", dir, cdffiles[i] ); ncid[i] = ncopen(filename,NC_NOWRITE); } /* * Get the lat/lon dimensions (they happen to be the * same for all files in this case) */ latid = ncdimid(ncid[0],"lat"); lonid = ncdimid(ncid[0],"lon"); ncdiminq(ncid[0],latid,(char *)0,&latlen); ncdiminq(ncid[0],lonid,(char *)0,&lonlen); len_dims[0] = latlen; len_dims[1] = lonlen; /* * Get the variable ids */ uid = ncvarid(ncid[0],"u"); vid = ncvarid(ncid[1],"v"); pid = ncvarid(ncid[2],"p"); tid = ncvarid(ncid[3],"t"); u5id = ncvarid(ncid[4],"u"); v5id = ncvarid(ncid[5],"v"); latid = ncvarid(ncid[0],"lat"); lonid = ncvarid(ncid[0],"lon"); /* * allocate space for arrays */ X = (float *)malloc(sizeof(float)*latlen*lonlen); Y = (float *)malloc(sizeof(float)*latlen*lonlen); lat = (float *)malloc(sizeof(float)*latlen); lon = (float *)malloc(sizeof(float)*lonlen); /* * Get lat/lon values (they are the same for all files) */ strt[0] = 0; cnt[0] = latlen; ncvarget(ncid[0],latid,(long const *)strt,(long const *)cnt,lat); cnt[0] = lonlen; ncvarget(ncid[0],lonid,(long const *)strt,(long const *)cnt,lon); /* * Get U and V data values */ get_2d_array(X,latlen,lonlen,ncid[0],uid,0); get_2d_array(Y,latlen,lonlen,ncid[1],vid,0); /* * Create a VectorField of the surface wind data */ NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNvfUDataArray,X,2,len_dims); NhlRLSetMDFloatArray(rlist,NhlNvfVDataArray,Y,2,len_dims); NhlRLSetFloat(rlist,NhlNvfXCStartV,lon[0]); NhlRLSetFloat(rlist,NhlNvfYCStartV,lat[0]); NhlRLSetFloat(rlist,NhlNvfXCEndV,lon[lonlen-1]); NhlRLSetFloat(rlist,NhlNvfYCEndV,lat[latlen-1]); NhlRLSetFloat(rlist,NhlNvfMissingUValueV,-9999.0); NhlCreate(&vfield,"VectorField",NhlvectorFieldClass,appid,rlist); /* * Create a VectorField of 500 millibar wind data * * Get U and V values */ get_2d_array(X,latlen,lonlen,ncid[4],u5id,0); get_2d_array(Y,latlen,lonlen,ncid[5],v5id,0); NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNvfUDataArray,X,2,len_dims); NhlRLSetMDFloatArray(rlist,NhlNvfVDataArray,Y,2,len_dims); NhlRLSetFloat(rlist,NhlNvfXCStartV,lon[0]); NhlRLSetFloat(rlist,NhlNvfYCStartV,lat[0]); NhlRLSetFloat(rlist,NhlNvfXCEndV,lon[lonlen-1]); NhlRLSetFloat(rlist,NhlNvfYCEndV,lat[latlen-1]); NhlRLSetFloat(rlist,NhlNvfMissingUValueV,-9999.0); NhlCreate(&vfield2,"VectorField",NhlvectorFieldClass,appid,rlist); /* * Create a ScalarField of surface pressure * * Get P data values */ get_2d_array(X,latlen,lonlen,ncid[2],pid,0); for( i = 0; i < latlen*lonlen; i++ ) { if( X[i] != -9999.0 ) { X[i] /= 100.; } } NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNsfDataArray,X,2,len_dims); NhlRLSetFloat(rlist,NhlNsfXCStartV,lon[0]); NhlRLSetFloat(rlist,NhlNsfYCStartV,lat[0]); NhlRLSetFloat(rlist,NhlNsfXCEndV,lon[lonlen-1]); NhlRLSetFloat(rlist,NhlNsfYCEndV,lat[latlen-1]); NhlRLSetFloat(rlist,NhlNsfMissingValueV,-9999.0); NhlCreate(&sfield,"ScalarField",NhlscalarFieldClass,appid,rlist); /* * Create a ScalarField of surface temperature * (convert from Kelvin to Farenheit) * * Get T data values */ get_2d_array(X,latlen,lonlen,ncid[3],tid,0); /* * Convert to Fahrenheit */ for( i = 0; i < latlen*lonlen; i++ ) { if( X[i] != -9999.0) { X[i] = (X[i] - 273.15) * 9.0/5.0 + 32.0; } } NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNsfDataArray,X,2,len_dims); NhlRLSetFloat(rlist,NhlNsfXCStartV,lon[0]); NhlRLSetFloat(rlist,NhlNsfYCStartV,lat[0]); NhlRLSetFloat(rlist,NhlNsfXCEndV,lon[lonlen-1]); NhlRLSetFloat(rlist,NhlNsfYCEndV,lat[latlen-1]); NhlRLSetFloat(rlist,NhlNsfMissingValueV,-9999.0); NhlCreate(&sfield2,"ScalarField2",NhlscalarFieldClass,appid,rlist); /* * Create a ContourPlot with surface temperature data */ NhlRLClear(rlist); NhlRLSetString(rlist,NhlNcnFillOn,"true"); NhlRLSetString(rlist,NhlNcnLinesOn,"false"); NhlRLSetString(rlist,NhlNcnFillDrawOrder,"predraw"); NhlRLSetInteger(rlist,NhlNcnScalarFieldData,sfield2); NhlCreate(&cnid,"contourplot",NhlcontourPlotClass,wid,rlist); /* * Create a VectorPlot with the surface wind and pressure data */ NhlRLClear(rlist); NhlRLSetString(rlist,NhlNvcUseScalarArray,"true"); NhlRLSetInteger(rlist,NhlNvcVectorFieldData,vfield); NhlRLSetInteger(rlist,NhlNvcScalarFieldData,sfield); NhlCreate(&vcid,"vectorplot",NhlvectorPlotClass,wid,rlist); /* * Create a StreamlinePlot with 500 mb wind data */ NhlRLClear(rlist); NhlRLSetString(rlist,NhlNpmTitleDisplayMode,"always"); NhlRLSetString(rlist,NhlNtiMainFuncCode,"~"); NhlRLSetInteger(rlist,NhlNstVectorFieldData,vfield2); NhlCreate(&stid,"streamlineplot",NhlstreamlinePlotClass,wid,rlist); /* * Create an annotation used to explain the streamline data */ NhlCreate(&txid,"streamlineplotanno",NhltextItemClass,wid,0); amid = NhlAddAnnotation(stid,txid); /* * Create a map object */ NhlRLClear(rlist); NhlRLSetString(rlist,NhlNvpUseSegments,"true"); NhlCreate(&mpid,"mapplot",NhlmapPlotClass,wid,rlist); /* * Overlay everything on the MapPlot. The last object overlaid will * appear on top */ NhlAddOverlay(mpid,cnid,-1); NhlAddOverlay(mpid,vcid,-1); NhlAddOverlay(mpid,stid,-1); /* * Variables for manipulating the title string */ tmid = ncdimid(ncid[1],"timestep"); ncdiminq(ncid[1],tmid,(char *)0,&timelen); tmid = ncvarid(ncid[1],"timestep"); timestep = (int *)malloc(sizeof(int)*timelen); strt[0] = 0; cnt[0] = timelen; ncvarget(ncid[1],tmid,(long const *)strt,(long const *)cnt,timestep); sprintf( hour, "00"); sprintf( day, "05"); stdmid = ncdimid(ncid[1],"timelen"); ncdiminq(ncid[1], stdmid, (char *)0, &stid_len ); tmid = ncvarid(ncid[1],"reftime"); rftime = (char *)malloc((stid_len+1)*sizeof(char)); strt[0] = 0; cnt[0] = stid_len; ncvarget(ncid[1],tmid,(long const *)strt,(long const *)cnt,rftime); for( i = 0; i <= TIMESTEPS-1; i++ ) { if (i != 17 && i != 36 && i != 37) { /* * Figure out the hour and day from the timestep, convert to strings * and build the title string */ d = timestep[i] / 24 + 5; h = timestep[i] % 24; if (h > 9) { sprintf( hour, "%d", h ); } else { sprintf( hour, "0%d", h ); } if (d > 9) { sprintf(day, "%d", d ); } else { sprintf(day, "0%d", d ); } /* * Set the new title string */ strcpy(mainstring, rftime); sprintf(&mainstring[8], "%2s %2s:00", day, hour); printf("%s\n",mainstring); NhlRLClear(rlist); NhlRLSetString(rlist,NhlNtiMainString,mainstring); NhlSetValues(stid,rlist); /* * Modify the data objects with data for the current time step * * Get U and V values */ get_2d_array(X,latlen,lonlen,ncid[0],uid,i); get_2d_array(Y,latlen,lonlen,ncid[1],vid,i); NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNvfUDataArray,X,2,len_dims); NhlRLSetMDFloatArray(rlist,NhlNvfVDataArray,Y,2,len_dims); NhlSetValues(vfield,rlist); /* * Get U and V values */ get_2d_array(X,latlen,lonlen,ncid[4],u5id,i); get_2d_array(Y,latlen,lonlen,ncid[5],v5id,i); NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNvfUDataArray,X,2,len_dims); NhlRLSetMDFloatArray(rlist,NhlNvfVDataArray,Y,2,len_dims); NhlSetValues(vfield2,rlist); /* * Get P values */ get_2d_array(X,latlen,lonlen,ncid[2],pid,i); for( j = 0; j < latlen*lonlen; j++ ) { if( X[j] != -9999.0 ) { X[j] /= 100.; } } NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNsfDataArray,X,2,len_dims); NhlSetValues(sfield,rlist); /* * Get T values */ get_2d_array(X,latlen,lonlen,ncid[3],tid,i); /* * Convert to Fahrenheit */ for( j = 0; j < latlen*lonlen; j++ ) { if( X[j] != -9999.0) { X[j] = (X[j] - 273.15) * 9.0/5.0 + 32.0; } } NhlRLClear(rlist); NhlRLSetMDFloatArray(rlist,NhlNsfDataArray,X,2,len_dims); NhlSetValues(sfield2,rlist); /* * Draw the plot */ NhlDraw(mpid); NhlFrame(wid); } } /* * Destroy the workstation object and exit. */ NhlDestroy(wid); NhlClose(); exit(0); } /* * function for reading in 3-d array from netCDF * file and converting it to a 2-d array. */ void get_2d_array( float *array, long latlen, long lonlen, int fid, int aid, long timestep ) { long start[3], count[3]; start[0] = timestep; start[1] = start[2] = 0; count[0] = 1; count[1] = latlen; count[2] = lonlen; ncvarget(fid,aid,(long const *)start,(long const *)count,array); return; }