PRO idl_hw05_acp ; ***************************** ; ******* Aaron Paget ******* ; ***************************** ; Purpose: This will extract out some wind from files and ; compute the magnitude of a subset. Then it is written into ; a new file which is smaller and easier to deal with. ; ; ; Project descriptions: ; Due Tuesday, next week, before class. ; Create a NetCDF file of the magnitude of the wind over North America. ; May-Sept, Long term mean values, up to 600mb, inclusive. ; lat = [15, 70], Lon [-135, -60] ; Sounds difficult, but its only tedious. ; Read data from uwnd, vwnd. Compute magnitude, save. ; Create new file, Define mode: set dimensions, set attributes ; In data mode, load data. Lat-lon-alt-time. (Jan=0, time array is ; [4,5,6,7,8]. Level=[1000, 925, 850, 700, 600] ; Use same names as CDC. Magnitude is wspd ; Note that you can have >1 file open at a time. ; Easy sequence of GET id_in, value to PUT, id_out, value ; ; Read data. File names ufn = '/net/b/klotter/uwnd.mon.ltm.nc' ; Set up the lats and lons LatMin = 15.; negative for South LatMax = 70. LonMin = -135. ; negative for West LonMax = -60. ; Choose a pressure level and month ; 0 = Jan, 1=Feb, etc... Pres = [1000., 925., 850., 700., 600.] ; milibar levels used Month = [4,5,6,7,8] ; May,June,July,Aug,Sept map_limit = [LatMin, LonMin, LatMax, LonMax]; for maps ;; Uwind data read lun_uwnd = NCDF_OPEN(ufn) ; note that it is assumed that the files start at Jan 1948 ; ; Retrieving attributes (ATTGET) from the NetCDF file: ; NCDF_ATTGET, lun_uwnd, "uwnd","add_offset", uw_add NCDF_ATTGET, lun_uwnd, "uwnd","scale_factor", uw_sf NCDF_ATTGET, lun_uwnd, "title", title, /GLOBAL NCDF_ATTGET, lun_uwnd, "history", history, /GLOBAL NCDF_ATTGET, lun_uwnd, "description", description, /GLOBAL NCDF_ATTGET, lun_uwnd, "platform", platform, /GLOBAL NCDF_ATTGET, lun_uwnd, "Conventions", Conventions, /GLOBAL ; Read data. Assume that the lon lat time and level are same for both U/V NCDF_VARGET, lun_uwnd, "lon", lon NCDF_VARGET, lun_uwnd, "lat", lat NCDF_VARGET, lun_uwnd, "time", time NCDF_VARGET, lun_uwnd, "level", level NCDF_VARGET, lun_uwnd, "uwnd", uwind ; assuming that SF=1, add_offset=0 NCDF_CLOSE, lun_uwnd ; closes the filenow that we read in the information ; selection parameters used for plot. limit to just those of interest here ; LONGITUDES. Complicated due to wrap-around for Prime Meridian lon_ind = [ WHERE(lon GE LonMin+360 AND lon LE LonMax+360), $ WHERE(lon GE LonMin AND lon LE LonMax) ] ; index of longitudes used check = WHERE(lon_ind GE 0, amt) lon_ind = lon_ind[check] lon_used = lon[lon_ind] lon_ind2 = WHERE(lon_used GT 180, amt) ; wrap around longitudes to plot over PM IF amt GT 0 THEN lon_used[lon_ind2] -= 360. ; value of lon's used n_lon = N_ELEMENTS(lon_used) ; get the information we want for our calculations lat_ind = WHERE(lat GE LatMin AND lat LE LatMax, n_lat) lat_used = lat[lat_ind] prs_ind = WHERE(Pres EQ level, n_prs) time_ind = Month n_tim=N_ELEMENTS(Month) ; set up the wind arrays uwind_used = FLTARR(n_lon, n_lat, n_prs, n_tim) FOR i=0,n_lon-1 DO BEGIN FOR j=0,n_lat-1 DO BEGIN FOR k=0,n_prs-1 DO BEGIN FOR L=0,n_tim-1 DO BEGIN uwind_used[i,j,k,L] = uwind[lon_ind[i], lat_ind[j], prs_ind[k], time_ind[L]] ENDFOR ENDFOR ENDFOR ENDFOR ; test print to verify it is reading correctly ;PRINT, 'ORIGINAL SIZE ', SIZE(uwind) ;PRINT, 'NEW SIZE ', SIZE(uwind_used) ; ***************Starting the V component.*********** ; Everything is the same as above for comment information please see ; the U component above for the V component and comments. ; ; Read data. File names vfn = '/net/b/klotter/vwnd.mon.ltm.nc' map_limit = [LatMin, LonMin, LatMax, LonMax]; for maps ; vwind data read lun_vwnd = NCDF_OPEN(vfn) ; note that it is assumed that the files start at Jan 1948 ; Read data. Assume that the lon lat time and level are same for both U/V NCDF_VARGET, lun_vwnd, "lon", lon NCDF_VARGET, lun_vwnd, "lat", lat NCDF_VARGET, lun_vwnd, "time", time NCDF_VARGET, lun_vwnd, "level", level NCDF_VARGET, lun_vwnd, "vwnd", vwind ; assuming that SF=1, add_offset=0 NCDF_CLOSE, lun_vwnd ; selection parameters used for plot. limit to just those of interest here ; LONGITUDES. Complicated due to wrap-around for Prime Meridian lon_ind = [ WHERE(lon GE LonMin+360 AND lon LE LonMax+360), $ WHERE(lon GE LonMin AND lon LE LonMax) ] ; index of longitudes used check = WHERE(lon_ind GE 0, amt) lon_ind = lon_ind[check] lon_used = lon[lon_ind] lon_ind2 = WHERE(lon_used GT 180, amt) ; wrap around longitudes to plot over PM IF amt GT 0 THEN lon_used[lon_ind2] -= 360. ; value of lon's used n_lon = N_ELEMENTS(lon_used) lat_ind = WHERE(lat GE LatMin AND lat LE LatMax, n_lat) lat_used = lat[lat_ind] prs_ind = WHERE(Pres EQ level, n_prs) prs_used=pres[prs_ind] time_ind = Month n_tim=N_ELEMENTS(Month) vwind_used = FLTARR(n_lon, n_lat, n_prs, n_tim) FOR i=0,n_lon-1 DO BEGIN FOR j=0,n_lat-1 DO BEGIN FOR k=0,n_prs-1 DO BEGIN FOR L=0,n_tim-1 DO BEGIN vwind_used[i,j,k,L] = vwind[lon_ind[i], lat_ind[j], prs_ind[k], time_ind[L]] ENDFOR ENDFOR ENDFOR ENDFOR ; test print ; PRINT, 'ORIGINAL SIZE ', SIZE(vwind) ; PRINT, 'NEW SIZE ', SIZE(vwind_used) ;PRINT, 'I got here' ; create the magnitude array wspd=SQRT(ABS(uwind_used^2.+vwind_used^2.)) ; PRINT, SIZE(wspd) ; CONTOUR, wspd(*,*,0,0), C_LABEL=[0,1,2,3,4,5,6] ; CONTOUR, uwind_used(*,*,0,0), C_LABEL=[0,1,2,3,4,5,6], COLOR='ff'x, /NOERASE ; CONTOUR, vwind_used(*,*,0,0), C_LABEL=[0,1,2,3,4,5,6], COLOR='ff00'x, /NOERASE ; ; ; ******************** Making the .NC file ************************ ; ; I give credit to Holly and Amanda for helping me with this. ; ; ; Generating the NetCDF file and over-writing any existing file: ; id = NCDF_CREATE('wspd_acp.nc', /CLOBBER) ; ; Generating default fill data into the file: ; NCDF_CONTROL, id, /FILL ; ; Creating the dimensions of the variables: ; timeid = NCDF_DIMDEF(id, 'r', 5) latid = NCDF_DIMDEF(id, 'x', 23) lonid = NCDF_DIMDEF(id, 'y', 31) presid = NCDF_DIMDEF(id, 'z', 5) ; ; Defining variables: ; tid = NCDF_VARDEF(id, 'time', [timeid], /SHORT) latid = NCDF_VARDEF(id, 'lat', [latid], /SHORT) lonid = NCDF_VARDEF(id, 'lon', [lonid], /SHORT) levid = NCDF_VARDEF(id, 'level', [presid], /SHORT) windid = NCDF_VARDEF(id, 'wspd', [lonid,latid,presid,timeid]) ; ; Inputting attributes into the file (gathered from initial NetCDF file): ; NCDF_ATTPUT, id, 'title', 'Monthly longterm mean windspeed calculated from u and v winds from the NCEP Reanalysis', /GLOBAL NCDF_ATTPUT, id, 'history', 'Created by Aaron Paget using NOAA-CIRES CDC data', /GLOBAL NCDF_ATTPUT, id, 'description', 'Windspeed for May through September', /GLOBAL NCDF_ATTPUT, id, "wspd","add_offset", uw_add NCDF_ATTPUT, id, "wspd","scale_factor", uw_sf NCDF_ATTPUT, id, tid, 'units', 'monthly' NCDF_ATTPUT, id, latid, 'units', 'degrees' NCDF_ATTPUT, id, latid, 'long_name', 'latitude' NCDF_ATTPUT, id, lonid, 'units', 'degrees' NCDF_ATTPUT, id, lonid, 'long_name', 'longitude' NCDF_ATTPUT, id, levid, 'units', 'millibars' NCDF_ATTPUT, id, levid,'long_name', 'pressure' NCDF_ATTPUT, id, windid, 'units', 'm/s' NCDF_ATTPUT, id, windid, 'long_name', 'Magnitude of wind speed' ; ; Take file out of define mode to data mode: ; NCDF_CONTROL, id, /ENDEF ; ; Writing data into the file: ; NCDF_VARPUT, id, tid, time_ind NCDF_VARPUT, id, latid, lat_used NCDF_VARPUT, id, lonid, lon_used NCDF_VARPUT, id, levid, prs_used NCDF_VARPUT, id, windid, wspd ; Read the data back out: ;NCDF_VARGET, id, tid, time ;NCDF_VARGET, id, latid, lat ;NCDF_VARGET, id, lonid, lon ;NCDF_VARGET, id, aid, level ;NCDF_VARGET, id, wid, wspd ;!P.CHARSIZE = 2.5 ;!X.TITLE = 'Location' ;!Y.TITLE = STRING(ytitle) ; Convert from bytes to strings. ;!Z.TITLE = STRING(ztitle) + '!C' + STRING(subtitle) ; ; Close NetCDF file: ; NCDF_CLOSE, id ; Close the NetCDF file. ; ; Hurray, it is over. ; END