;---------------------------------------------------------------------- ; shapefiles_11.ncl ; ; Concepts illustrated: ; - Reading shapefiles ; - Creating a mask variable using outlines from a shapefile ; - Attaching markers to a map ; - Attaching polylines to a map plot ; - Masking a data array based on a geographical area obtained from a shapefile ;---------------------------------------------------------------------- ; This script creates a new data mask based on the outline of the ; United States. It then draws two plots: the original data, ; and the data with the USA mask. ; ; The "gz_2010_us_020_00_5m.shp" shapefile is from ; http://www.census.gov/geo/www/cob/ ; ; The "nationalp010g.shp" shapefile is from: ; http://www.nationalatlas.gov/atlasftp-1m.html#nationp ; ; Try both of these files to see if they work for your purposes. ; The national file is larger, and hence takes longer to process ; (206 wall clock seconds on a Mac versus 17 wall clock seconds.) ; ; Once you create the mask, you can set CREATE_MASK to False so ; it doesn't get created again when you run this script. ; ;---------------------------------------------------------------------- load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" load "shapefile_mask_data.ncl" ;---------------------------------------------------------------------- ; This is the main code ;---------------------------------------------------------------------- begin ;pltDir = ("/Amal/") ; graphic directory pltDir = ("/home/netapp-clima/scratch/asalhi/mask/") ; graphic directory ; ; If you already have the mask NetCDF file, set this to False. ; Creating the mask can be slow! ; CREATE_MASK = True ;---Whether to draw lat/lon points and shapefile outlines ADD_LATLON_POINTS = False ADD_SHAPEFILE_OUTLINES = True ;---Name of shapefile containing USA outlines shp_fname = "GSAs_simplified.shp" ;---Name of file to write mask to or to read mask from. mask_fname = "created_mask.nc" ;---Rough area we are interested in. Everything outside this will be masked. minlon = -08.0 maxlon = 45.0 minlat = 29.00 maxlat = 48.0 ; minlon = -10.0 ; maxlon = 40.0 ; minlat = 30.00 ; maxlat = 50.0 ;---Read in zonal winds a = addfile("bbio.nc","r") ;b = addfile("bbio-ts.nc","w") dims = getfiledimsizes(a) print(dims) nt = dims(3) ; print(nt) ; n = 0 ; do while (n.le.nt) ; f = a->dox(n,0,:,:) f = a->dox(0,0,:,:) if(CREATE_MASK) then print("Creating the mask file...") ;---Create a new mask using a shapefile of USA pdims = dimsizes(f) opt = True opt@return_mask = True ; opt@debug = True opt@minlon = minlon ; Makes the shapefile masking opt@maxlon = maxlon ; go faster. opt@minlat = minlat opt@maxlat = maxlat opt@shape_var = "SECT_COD" opt@shape_names = ("GSA21") oxygen = shapefile_mask_data(f,shp_fname,opt) copy_VarCoords(f,oxygen) printVarSummary(oxygen) printMinMax(oxygen,0) ;---Write new mask to file system("rm -f " + mask_fname) fout = addfile(mask_fname,"c") fout->oxygen_mask = oxygen else print("Reading mask off file.") ;---Read the new mask from the NetCDF file fmask = addfile(mask_fname,"r") oxygen = fmask->oxygen_mask end if ;---Create masked data array pmask = where(oxygen.eq.1,f,f@_FillValue) copy_VarMeta(f,pmask) ;---Start the graphics. plot = new(1,graphic) wks = gsn_open_wks("eps","oxygen_aa") ; wks = gsn_open_wks("eps","oxygen"+n) gsn_define_colormap(wks,"CBR_coldhot") gsn_reverse_colormap(wks) res = True ; plot mods desired res@gsnPaperOrientation = "portrait" ; landscape, portrait; Change orientation res@gsnMaximize = True res@gsnDraw = False res@gsnFrame = False res@gsnLeftString = "" ; don't draw left string; ne pas afficher le nom de la variable sur la figure res@gsnRightString = "" ; don't draw right string;ne pas afficher le nom de la variable sur la figure res@mpOutlineOn = False res@cnFillOn = True ; turn on color fill res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour line labels res@lbLabelBarOn = False ; will turn on in panel res@mpDataBaseVersion = "MediumRes" res@mpMinLonF = minlon res@mpMaxLonF = maxlon res@mpMinLatF = minlat res@mpMaxLatF = maxlat res@mpCenterLonF = (minlon+maxlon)/2. res@gsnAddCyclic = False ; Regional data ;---Be sure to use same levels for both plots ;---Create (but don't draw) both plots ;res@tiMainString = "Original precip" ;plot_orig = gsn_csm_contour_map(wks,f,res) ;res@tiMainString = "Dissolved_Molecular_Oxygen" res@gsnRightString = "mmol/m3" res@gsnLeftString = "Dissolved_Molecular_Oxygen" res@cnLevelSelectionMode = "ManualLevels" ; manually set the contour levels with the following 3 resources res@cnMinLevelValF = 0. ; set the minimum contour level res@cnMaxLevelValF = 300. ; set the maximum contour level res@cnLevelSpacingF = 10 ;res@tiMainString = "Masked precip from " + shp_fname plot_mask = gsn_csm_contour_map(wks,pmask,res) if(ADD_LATLON_POINTS) then ;---Set up a resource list to attach the grid points as filled dots ;mkres1 = True ;mkres1@gsnCoordsAttach = True ;gsn_coordinates(wks,plot_orig,f,mkres1) mkres2 = True mkres2@gsnCoordsAttach = True mkres2@gsnCoordsNonMissingColor = "black" mkres2@gsnCoordsMissingColor = "red" gsn_coordinates(wks,plot_mask,pmask,mkres2) end if if(ADD_SHAPEFILE_OUTLINES) then lnres = True ;poly_orig = gsn_add_shapefile_polylines(wks,plot_orig,shp_fname,lnres) poly_mask = gsn_add_shapefile_polylines(wks,plot_mask,shp_fname,lnres) end if ;---Draw both plot in a panel. pres = True pres@gsnMaximize = True pres@gsnPanelLabelBar = True pres@pmLabelBarWidthF = 0.8 ;pres@gsnPanelFigureStrings= (/"a)","b)"/) ; add strings to panel pres@gsnPanelFigureStringsFontHeightF = 0.02 ; decrease the size of figure strings pres@amJust = "TopLeft" ; Where to put these strings pres@gsnPanelYWhiteSpacePercent = 2 ; Adds the white space to the panel plot pres@gsnPanelXWhiteSpacePercent = 2 ; Adds the white space to the panel plot pres@txFontHeightF = 0.025 ; smaller text font size ;pres@txString = "pr_masked" ; add main title for the panel ; Changing the font seize pres@tmXBLabelFontHeightF = 0.018 ; resize tick labels for X axis pres@tmYLLabelFontHeightF = 0.018 ; resize tick labels for Y axis ; change label spacing to avoid overlap pres@lbLabelStride = 2 ; and write every other label pres@lbLabelFontHeightF = 0.018 ; change the size of the label bar pres@lbBoxMinorExtentF = 0.15 ;-- decrease the height of the labelbar pres@IbLabelFontHeightF = 0.01 ; Decreases the label sizes of the common labelbar pres@gsnStringFontHeightF = 0.02 ; Increase the title font size (long_name) gsn_panel(wks,(/plot_mask/),(/1,1/),pres) ;----create ncfile ntim = nt ; get dimension sizes nlev = dims(2) nlat = dims(1) nlon = dims(0) ; nvar =1 time=a->time lat=a->latitude lon=a->longitude lev=a->depth diro = "./" ; Output directory filo = "testaa.nc" ; filo = fname ; Output file system("/bin/rm -f " + diro + filo) ; remove if exists foutput = addfile (diro + filo, "c") ; open output file setfileoption(foutput,"DefineMode",True) fAtt = True ; assign file attributes fAtt@title = "NCL Efficient Approach to netCDF Creation" fAtt@source_file = "bbio.nc" fAtt@Conventions = "None" fAtt@creation_date = systemfunc ("date") fileattdef( foutput, fAtt ) ; copy file attributes dimNames = (/"time","depth", "latitude", "longitude"/) dimSizes = (/ -1 , nlev, nlat, nlon /) dimUnlim = (/ True , False, False, False/) filedimdef(foutput,dimNames,dimSizes,dimUnlim) oxy=pmask print("time="+time) ;pmask=new((/nt,nlev,nlat,nlon/), float) ; oxy@time = time ; Name the dimensions ; oxy@lev = lev pmask@lat = lat pmask@lon = lon filevardef(foutput, "time" ,typeof(time),getvardims(time)) filevardef(foutput, "lev" ,typeof(lev),getvardims(lev) ) filevardef(foutput, "lat" ,typeof(lat),getvardims(lat)) filevardef(foutput, "lon" ,typeof(lon),getvardims(lon)) filevardef(foutput, "pmask" ,typeof(pmask) ,getvardims(pmask)) ; filevardef(fout, "PS" ,typeof(PS) ,getvardims(PS)) ; filevardef(fout, "TOPOG",typeof(ORO),getvardims(ORO)) ; variable name on the file ; different from name on script filevarattdef(foutput,"pmask",pmask) ; copy T attributes filevarattdef(foutput,"time" ,time) ; copy time attributes filevarattdef(foutput,"lev" ,lev) ; copy lev attributes filevarattdef(foutput,"lat" ,lat) ; copy lat attributes filevarattdef(foutput,"lon" ,lon) ; copy lon attributes ; filevarattdef(fout,"PS" ,PS) ; copy PS attributes ; filevarattdef(fout,"TOPOG",ORO) ; copy TOPOG attributes setfileoption(foutput,"DefineMode",False) foutput->time = (/time/) foutput->lev = (/lev/) foutput->lat = (/lat/) foutput->lon = (/lon/) foutput->pmask = (/pmask/) ; fout->Z = (/PS/) ; fout->TOPOG = (/ORO/) end