#!/usr/bin/env python # ============================================================================ # MEDESS_FORC_ORIG.PY # # This software read NetCDF Medslik-II outputs and plots: # # 1. Oil concentration fields on a 150x150 m regular grid # # 2. Ocean currents velocities fields on the original ocean currents grid # # 3. Stokes' Drift velocities fields on the original waves grid # #------------------------------------------------------------------------------ import sys import numpy as np import os import Nio import Ngl import mpl_toolkits.basemap.pyproj as pyproj import shapefile dir = sys.argv[1] file = sys.argv[2] name = dir + file f_nc = Nio.open_file(name,"r") cst_file = dir + "runcoast.map" # ---------------------------------------------------------------- # TIME variable # ---------------------------------------------------------------- time = f_nc.variables["time"] # ---------------------------------------------------------------- # OIL CONCENTRATION VARIABLE # ---------------------------------------------------------------- oil = f_nc.variables["surface"] #print "oil shape = ", oil.shape #print "oil dimensions = ", oil.dimensions #oilAttsAndVals = oil.attributes #print "Dictionary of Attibutes and Values for oil: ", oilAttsAndVals if hasattr(oil,"_FillValue"): NaN_oil = oil._FillValue print "oil FillValue = ", NaN_oil if hasattr(oil,"scale_factor"): sc_f_oil = oil.scale_factor print "oil scale_factor = ", sc_f_oil else: sc_f_oil = 1. if hasattr(oil,"add_offset"): ad_o_oil = oil.add_offset print "conc add offset = ", ad_o_oil else: ad_o_oil = 0. C1 = (ad_o_oil + (oil[1:len(time[:]),:,:] * sc_f_oil)) #print C1.dtype #print C1.shape Cmin = np.amin(C1) print "Cmin: ", Cmin Cmax = np.amax(C1) print "Cmax: ", Cmax C = (ad_o_oil + (oil[:,:,:] * sc_f_oil)) # ---------------------------------------------------------------- # OIL ON COAST CONCENTRATION VARIABLE # ---------------------------------------------------------------- cst_seg = open(cst_file).readlines() seg = np.zeros((len(cst_seg),6)) if len(cst_seg) != 0: for i in range(0,len(cst_seg)): if float(cst_seg[i].split()[0]) != 999999: for j in range(0,6): seg[i,j] = float(cst_seg[i].split()[j]) varNames = f_nc.variables.keys() if varNames[9] == "coast": cst = f_nc.variables["coast"] print "cst shape = ", cst.shape print "cst dimensions = ", cst.dimensions cstAttsAndVals = cst.attributes print "Dictionary of Attibutes and Values for oil on coast: ", cstAttsAndVals if hasattr(cst,"_FillValue"): NaN_cst = cst._FillValue print "cst FillValue = ", NaN_cst if hasattr(cst,"scale_factor"): sc_f_cst = cst.scale_factor print "cst scale_factor = ", sc_f_cst else: sc_f_cst = 1. if hasattr(cst,"add_offset"): ad_o_cst = cst.add_offset print "conc add offset = ", ad_o_cst else: ad_o_cst = 0. CST1 = (ad_o_cst + (cst[1:len(time[:]),:] * sc_f_cst)) print CST1.dtype print CST1.shape COASTmin = np.amin(CST1) print "COASTmin: ", COASTmin COASTmax = np.amax(CST1) print "COASTmax: ", COASTmax CST = (ad_o_cst + (cst[:,:] * sc_f_cst)) else: CST = np.zeros((len(time),1)) # --------------------------------------------------------------- # CURRENTS U and V variables # --------------------------------------------------------------- uu = f_nc.variables["u_cur"] vv = f_nc.variables["v_cur"] #print "u shape = ", uu.shape #print "u dimensions = ", uu.dimensions #print "v shape = ", vv.shape #print "v dimensions = ", vv.dimensions #uAttsAndVals = uu.attributes #vAttsAndVals = vv.attributes #print "Dictionary of Attibutes and Values for u: ", uAttsAndVals #print "Dictionary of Attibutes and Values for v: ", vAttsAndVals if hasattr(uu,"_FillValue"): NaN_c = uu._FillValue print "u FillValue = ", NaN_c if hasattr(uu,"scale_factor"): sc_f_c = uu.scale_factor print "conc scale_factor = ", sc_f_c else: sc_f_c = 1. if hasattr(uu,"add_offset"): ad_o_c = uu.add_offset print "u add offset = ", ad_o_c else: ad_o_c = 0. u = (ad_o_c + (uu[:,:,:] * sc_f_c)) v = (ad_o_c + (vv[:,:,:] * sc_f_c)) print u.dtype, v.dtype print u.shape, v.shape umin = np.amin(u) print "umin: ", umin umax = np.amax(u) print "umax: ", umax vmin = np.amin(v) print "vmin: ", vmin vmax = np.amax(u) print "vmax: ", vmax # --------------------------------------------------------------- # WAVE U and V variables # --------------------------------------------------------------- wu = f_nc.variables["u_wav"] wv = f_nc.variables["v_wav"] #print "u shape = ", wu.shape #print "u dimensions = ", wu.dimensions #print "v shape = ", wv.shape #print "v dimensions = ", wv.dimensions #uwAttsAndVals = wu.attributes #vwAttsAndVals = wv.attributes #print "Dictionary of Attibutes and Values for u: ", uwAttsAndVals #print "Dictionary of Attibutes and Values for v: ", vwAttsAndVals if hasattr(wu,"_FillValue"): NaN_w = wu._FillValue print "u FillValue = ", NaN_w if hasattr(wu,"scale_factor"): sc_f_w = wu.scale_factor print "conc scale_factor = ", sc_f_w else: sc_f_w = 1. if hasattr(wu,"add_offset"): ad_o_w = wu.add_offset print "u add offset = ", ad_o_w else: ad_o_w = 0. uw = (ad_o_w + (wu[:,:,:] * sc_f_w)) vw = (ad_o_w + (wv[:,:,:] * sc_f_w)) print uw.dtype, vw.dtype print uw.shape, vw.shape uwmin = np.amin(uw) print "uwmin: ", uwmin uwmax = np.amax(uw) print "uwmax: ", uwmax vwmin = np.amin(vw) print "vwmin: ", vwmin vwmax = np.amax(vw) print "vwmax: ", vwmax # ------------------------------------------------------------------------------- # OIL LON and LAT variables # ------------------------------------------------------------------------------- lon = f_nc.variables["lon_oil"] lat = f_nc.variables["lat_oil"] print "lon shape = ", lon.shape print "lon dimensions = ", lon.dimensions print "lat shape = ", lat.shape print "lat dimensions = ", lat.dimensions nlon = len(lon[:]) nlat = len(lat[:]) # ------------------------------------------------------------------------------- # CURRENTS LON and LAT variables # ------------------------------------------------------------------------------- lon_cur = f_nc.variables["lon_cur"] lat_cur = f_nc.variables["lat_cur"] print "lon_cur shape = ", lon_cur.shape print "lon_cur dimensions = ", lon_cur.dimensions #print lon_cur.get_value() print "lat_cur shape = ", lat_cur.shape print "lat_cur dimensions = ", lat_cur.dimensions #print lat_cur.get_value() nlon_cur = len(lon_cur[:]) nlat_cur = len(lat_cur[:]) # ------------------------------------------------------------------------------- # WAVES LON and LAT variables # ------------------------------------------------------------------------------- lon_wav = f_nc.variables["lon_wav"] lat_wav = f_nc.variables["lat_wav"] print "lon_wav shape = ", lon_wav.shape print "lon_wav dimensions = ", lon_wav.dimensions #print lon_wav.get_value() print "lat_wav shape = ", lat_wav.shape print "lat_wav dimensions = ", lat_wav.dimensions #print lat_wav.get_value() nlon_wav = len(lon_wav[:]) nlat_wav = len(lat_wav[:]) #--------------------------------------------------------------------------------- # EMODNET NATURA2000 SHAPEFILE #--------------------------------------------------------------------------------- file = "/home/d13g0/MEDSLIKII_plot/python_code/python_shapefile/Natura2000_end2013_rev1_Shapefile/Natura2000_end2013_rev1_med.shp" etrs89=pyproj.Proj("+init=EPSG:3035") # ETRS89 / ETRS-LAEA, single Coordinate Reference System for all Europe sf = Nio.open_file(file,"r") # Read variables geom = sf.variables["geometry"][:] # geom is a numpy matrix segm = sf.variables["segments"][:] # segm is a numpy matrix lonkm = sf.variables["x"][:] print type(lonkm) latkm = sf.variables["y"][:] LON = np.zeros((len(lonkm))) LAT = np.zeros((len(latkm))) print "DONE READ VARIABLES" LON, LAT = etrs89(lonkm,latkm, inverse=True) print "DONE CONVERSION OF COORDINATE REFERENCE SYSTEM" # Read variables attributes geomDims = geom.shape print "geomDims: ", geomDims segsDims = segm.shape print "segsDims: ", segsDims numFeatures = geomDims[0] # Read global attributes geom_segIndex = sf.geom_segIndex print "geom_segIndex: ", geom_segIndex geom_numSegs = sf.geom_numSegs print "geom_numSegs: ", geom_numSegs segs_xyzIndex = sf.segs_xyzIndex print "segs_xyzIndex: ", segs_xyzIndex segs_numPnts = sf.segs_numPnts print "segs_numPnts: ", segs_numPnts print "START TO PLOT" # -------------------------------------------------------------------------------- # LOOP on time to plot scalar and vector fields # -------------------------------------------------------------------------------- for t in range(0, len(time[:])): print "TIMESTEP: ", t C1 = C[t,:,:] CST1 = CST[t,:] ut = u[t,:,:] vt = v[t,:,:] u_log = ut == 9999. v_log = vt == 9999. ut[u_log] = 0. vt[v_log] = 0. uwt = uw[t,:,:] vwt = vw[t,:,:] uw_log = uwt == 9999. vw_log = vwt == 9999. uwt[uw_log] = 0. vwt[vw_log] = 0. numin = np.amin(ut) numax = np.amax(ut) nvmin = np.amin(vt) nvmax = np.amax(vt) nuwmin = np.amin(uwt) nuwmax = np.amax(uwt) nvwmin = np.amin(vwt) nvwmax = np.amax(vwt) # Plot ------------------------- # 1. Open the WorKStation out_name = "medslikII_" + str(t) rlist = Ngl.Resources() rlist.wkColorMap = "WhBlGrYeRe" #cmap wks_type = "png" wks = Ngl.open_wks(wks_type,out_name,rlist) cmap = Ngl.retrieve_colormap(wks) # print cmap # 2. Set up resource lists will apply to vector, line contour, and # filled contour plots. mpres = Ngl.Resources() # map resources c_vcres = Ngl.Resources() # ocean currents vector resources w_vcres = Ngl.Resources() # Stokes' Drift vector resources cfres = Ngl.Resources() # contour line/fill resources cstres = Ngl.Resources() # Coastal segmments polyline resources natres = Ngl.Resources() # Natura2000 segmments polyline resources # 3. Turn off nglDraw and nglFrame because we don't want to draw all # these plots until they are all overlaid on the map plot. mpres.nglDraw = False mpres.nglFrame = False c_vcres.nglDraw = False c_vcres.nglFrame = False w_vcres.nglDraw = False w_vcres.nglFrame = False cfres.nglDraw = False cfres.nglFrame = False cstres.nglDraw = False cstres.nglFrame = False natres.nglDraw = False natres.nglFrame = False # 4. Set up coordinates of X and Y axes for all plots. This # is necessary in order for the Ngl.overlay calls to work # later. cfres.sfXCStartV = float(lon[0]) # Define X/Y axes range cfres.sfXCEndV = float(lon[nlon-1]) # for contour plot. cfres.sfYCStartV = float(lat[0]) cfres.sfYCEndV = float(lat[nlat-1]) c_vcres.vfXCStartV = float(lon_cur[0]) # Define X/Y axes range c_vcres.vfXCEndV = float(lon_cur[nlon_cur-1]) # for vector plot. c_vcres.vfYCStartV = float(lat_cur[0]) c_vcres.vfYCEndV = float(lat_cur[nlat_cur-1]) w_vcres.vfXCStartV = float(lon_wav[0]) # Define X/Y axes range w_vcres.vfXCEndV = float(lon_wav[nlon_wav-1]) # for vector plot. w_vcres.vfYCStartV = float(lat_wav[0]) w_vcres.vfYCEndV = float(lat_wav[nlat_wav-1]) # 5. Set additional resources needed. # CURRENTS VECTOR FIELD SETUP------------------------------------------------------------ # vcres.vcMinFracLengthF = 0.33 # Increase length of # vcres.vcMinMagnitudeF = 0.001 # vectors. c_vcres.vcGlyphStyle = "LineArrow" c_vcres.vcMonoLineArrowColor = True c_vcres.vcLineArrowColor = "black" c_vcres.vcRefLengthF = 0.05 c_vcres.vcVectorDrawOrder = "PreDraw" c_vcres.vcRefMagnitudeF = 0.2 c_vcres.vcPositionMode = "ArrowTail" c_vcres.vcLineArrowThicknessF = 2. # Double the thickness. # c_vcres.vcRefAnnoOrthogonalPosF = -0.179 # Move reference annotation up c_vcres.vcRefAnnoParallelPosF = 0.19 # and over to left. c_vcres.vcRefAnnoString1 = "0.2 m/s" c_vcres.vcRefAnnoString2 = "Ocean Currents" # WAVES VECTOR FIELD SETUP------------------------------------------------------------ # vcres.vcMinFracLengthF = 0.33 # Increase length of # vcres.vcMinMagnitudeF = 0.001 # vectors. w_vcres.vcGlyphStyle = "LineArrow" w_vcres.vcMonoLineArrowColor = True w_vcres.vcLineArrowColor = "green" w_vcres.vcRefLengthF = 0.05 w_vcres.vcVectorDrawOrder = "PreDraw" w_vcres.vcRefMagnitudeF = 0.1 w_vcres.vcPositionMode = "ArrowTail" w_vcres.vcLineArrowThicknessF = 2. # Double the thickness. # vcres.vcRefAnnoOrthogonalPosF = -0.179 # Move reference annotation up # vcres.vcRefAnnoParallelPosF = 0.19 # and over to left. w_vcres.vcRefAnnoString1 = "0.1 m/s" w_vcres.vcRefAnnoString2 = "Wave-induced Currents" # SCALAR FIELD SETUP cfres.cnFillOn = True # Turn on contour fill. cfres.cnLinesOn = False # Turn on contour lines. cfres.cnLineLabelsOn = False # Turn off contour line labels. cfres.lbOrientation = "Horizontal" # horizontal labelbar cfres.lbLabelFontHeightF = 0.012 # Decrease font size. cfres.pmLabelBarOrthogonalPosF = -0.05 # Move labelbar up. cfres.cnLineDashPatterns = 3 # dashed contour lines cfres.cnLineThicknessF = 3. # triple thick contour lines cfres.cnInfoLabelOrthogonalPosF = -0.15 # Move info label up. cfres.cnFillMode = "RasterFill" # or "AreaFill" or "CellFill" or "RasterFill" cfres.cnRasterSmoothingOn = True cfres.cnLineDrawOrder = "PreDraw" cfres.cnFillDrawOrder = "PreDraw" cfres.cnLevelSelectionMode = "ExplicitLevels" # Define your own contour levels cfres.cnLevels = np.arange(0.001,Cmax,0.001) # cfres.cnLevels = np.arange(0.1,Cmax,0.1) # cfres.cnLevels = np.arange(0.001,0.015,0.001) cfres.lbOrientation = "Horizontal" cfres.lbTitleString = "Surface oil concentration, Kg/m^2" cfres.lbTitleFontHeightF = 0.012 cfres.lbLabelFontHeightF = 0.008 cfres.lbTitleOffsetF = -0.27 cfres.lbBoxMinorExtentF = 0.15 cfres.pmLabelBarOrthogonalPosF = -0.01 # POLYLINE SETUP cstres.gsLineThicknessF = 8.0 cstres.gsLineColor = "black" natres.gsLineColor = "red" # 6. Draw VECTOR and SCALAR FIELD cf = Ngl.contour(wks,C1,cfres) vccur = Ngl.vector(wks,ut,vt,c_vcres) vcwav = Ngl.vector(wks,uwt,vwt,w_vcres) # 7. MAP SETUP xs = Ngl.get_float(cf.sffield,"sfXCActualStartF") - 0.2 xe = Ngl.get_float(cf.sffield,"sfXCActualEndF") + 0.2 ys = Ngl.get_float(cf.sffield,"sfYCActualStartF") - 0.2 ye = Ngl.get_float(cf.sffield,"sfYCActualEndF") + 0.2 # print xs, xe # print ys, ye # xs = 7.50 # xe = 8.50 # ys = 39.00 # ye = 40.00 mpres.mpProjection = "CylindricalEquidistant" mpres.mpDataBaseVersion = "HighRes" mpres.mpLimitMode = "LatLon" mpres.mpMinLonF = xs mpres.mpMaxLonF = xe mpres.mpMinLatF = ys mpres.mpMaxLatF = ye mpres.mpPerimOn = True # Turn on map perimeter. mpres.mpGridAndLimbOn = True mpres.mpPerimDrawOrder = "PostDraw" mpres.mpFillDrawOrder = "PostDraw" mpres.tfPolyDrawOrder = "PostDraw" mpres.mpOutlineBoundarySets = "GeophysicalAndUSStates" mpres.mpFillOn = True mpres.mpFillColors = ["background","transparent","LightGray","LightGray"] # Fill land and inland water # and leave ocean transparent # mpres.mpFillColors = ["background","transparent","transparent","transparent"] mpres.pmTitleDisplayMode = "Always" # Turn on map title. mpres.pmTickMarkDisplayMode = "Always" mpres.pmLabelBarOrthogonalPosF = -0.05 # 8. TITLE resources mpres.tiMainString = "MEDSLIK-II SIM 2" mpres.tiMainFontHeightF = 0.020 mpres.tiMainFont = "Helvetica-bold" mpres.tiMainOffsetYF = 0.025 # 9. Plot the map mp = Ngl.map(wks,mpres) # Draw a plot of map # 8. Overlay everything on the map plot. Ngl.overlay(mp,cf) Ngl.overlay(mp,vccur) Ngl.overlay(mp,vcwav) # 10. Add Natura2000 polylines to map for i in range(0,numFeatures): startSegment = geom[i, geom_segIndex] numSegments = geom[i, geom_numSegs] for seg in range(startSegment, startSegment+numSegments): startPT = segm[seg, segs_xyzIndex] endPT = startPT + segm[seg, segs_numPnts]# - 1 Ngl.polyline(wks, mp, LON[startPT:endPT], LAT[startPT:endPT], natres) # 9. Add the Coast segments if present for i in range(0,len(CST1)): if CST1[i] != 0.: x1 = seg[i,0] x2 = seg[i,2] y1 = seg[i,1] y2 = seg[i,3] Ngl.add_polyline(wks,mp,[x1,x2],[y1,y2],cstres) # 11. Draw the map plot, which now contains the vectors and # filled/line contours. Ngl.maximize_plot(wks,mp) # Maximize size of plot in frame. Ngl.draw(mp) # # Draw text strings after the plot is drawn to make sure plot # gets maximized properly. # txres = Ngl.Resources() txres.txFontHeightF = 0.015 count_date = t + 6 if (count_date <= 24): string1 = "17-05-2014" else: string1 = "18-05-2014" if (t < 10) : string2 = "+ 000" + str(t) else: string2 = "+ 00" + str(t) string = string1 + " 05:38 UTC " + string2 + " hours " Ngl.text_ndc(wks,string, 0.335, 0.920, txres) # Ngl.text_ndc(wks,string2, 0.85, 0.915, txres) # Ngl.text_ndc(wks,"Max Elevation: 4322", 0.85, 0.775, txres) Ngl.frame(wks) Ngl.destroy(wks) Ngl.end()