Re: Trying to get correct map over rotated lat-lon data

Dear Dave, thank you so much for your further suggestions, time and effort. These latest comments have indeed made improvements and the script is now working in a way that I would consider to be more than acceptable. Once again the NCL-talk team has gone beyond the call of duty to provide an excellent service. You have saved me many headaches and your efforts are very much appreciated. I will attach and paste below a cleaned up version of my working script in case others have this problem in the future. Yours in gratitude, Dave. -----------SCRIPT--------------------------------- ; Basic script for plotting model output 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/wrf/WRFUserARW.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin a = addfile("wrfout_d01_2009-12-01_00:00:00.nc","r") wks = gsn_open_wks("x11","WRF_RotLatLon") gsn_define_colormap(wks,"precip3_16lev") ;--------------------------------------------------------------- ; Set basic resources mpres = True mpres_at_tiMainString = "Main Title"        ; Give a main title ;--------------------------------------------------------------- ; Get variables ;---1. WRF------------------ ; ;psfc          = a->PSFC(0,:,:)        ; switch between psfc or landmask psfc          = a->LANDMASK(0,:,:)  ; to test fit of geo borders psfc_size     = dimsizes(psfc(:,:)) print("dimsizes(psfc) = " + psfc_size) ;---2. 2D lat/lon arrays---- lat2d  = a->XLAT(0,:,:) lon2d  = a->XLONG(0,:,:) ; Assign lat2d & lon2d attributes to variable psfc_at_lat2d   = True psfc_at_lon2d   = True ; could alternatively set mpres_at_gsnAddCyclic=False ; Get min/max extent of domain grid points max_y  = dimsizes(lat2d(:,0)) - 1 max_x  = dimsizes(lon2d(0,:)) - 1 print("max_y = " + max_y) print("max_x = " + max_x) ; Set some map options     mpres_at_tfDoNDCOverlay= True ; Important!!     mpres_at_mpGridAndLimbOn= True     mpres_at_mpGridSpacingF = 10     mpres_at_mpCenterLonF   = -a_at_STAND_LON        mpres_at_mpCenterLatF   = 90 - a_at_POLE_LAT     mpres_at_mpLimitMode = "corners"     mpres_at_mpLeftCornerLonF = lon2d(0,0)             mpres_at_mpRightCornerLonF= lon2d(max_y,max_x)     mpres_at_mpLeftCornerLatF = lat2d(0,0)            mpres_at_mpRightCornerLatF= lat2d(max_y,max_x)        mpres_at_mpProjection   = "CylindricalEquidistant"             mpres_at_cnFillOn    = True     mpres_at_cnLinesOn    = False     mpres_at_gsnMaximize = True         mpres_at_mpOutlineOn = True     mpres_at_mpGeophysicalLineColor = "Black"     mpres_at_mpGridLineColor = "Black"     mpres_at_mpLimbLineColor = "Black"     mpres_at_mpNationalLineColor = "Black"     mpres_at_mpUSStateLineColor = "Black"     mpres_at_mpPerimLineColor = "Black"     mpres_at_mpGeophysicalLineThicknessF = 1.0     mpres_at_mpGridLineThicknessF = 1.0     mpres_at_mpLimbLineThicknessF = 1.0     mpres_at_mpNationalLineThicknessF = 1.0     mpres_at_mpUSStateLineThicknessF = 1.0 ;  Namelist map projection information: ;  map_proj : "lat-lon" (cylindrical equidistant) ;  ref_lat  : (geographical lat of grid point (ref_x,ref_y) ;  ref_lon  : (geographical lon of grid point (ref_x,ref_y) ;  ref_x    : (x grid point) ;  ref_y    : (y grid point) ;  pole_lat : (lat of geog N hem pole on comp grid) ;  pole_lon : (lon of geog N hem pole on comp grid) ;  stand_lon: (rotation about geographic poles) ; MAKE PLOT  plot = gsn_csm_contour_map_ce(wks,psfc,mpres)   delete(mpres) end ________________________________ From: David Brown <dbrown_at_ucar.edu> To: David Jones <jonesd647_at_yahoo.ca> Cc: ncl-talk Talk <ncl-talk_at_ucar.edu> Sent: Friday, November 30, 2012 8:42:47 PM Subject: Re: Trying to get correct map over rotated lat-lon data Hi David, The issue with sfXArray and sfYArray is that, by default, the gsn_csm routines set the resource gsnAddCyclic to True. This appends a copy of the first column of the data to the end of the array resulting in the X dimension of the data array increasing by one element. I have to note here that there are 2 ways to supply 2d coordinate arrays to the gsn routines. You can 1) attach the special attributes lon2d and lat2d to the data array, or 2) you can set sfXArray and sfYArray. In the first case, the gsn routines set sfXArray and sfYArray for you using the values supplied as attributes, and also make sure to add the extra column to both the coordinate arrays. Unfortunately they do not adjust the sfXArrary and sfYArray values if you set them directly.  The purpose of the additional column is to eliminate the gap that can appear in plots of global grids at the cyclic longitude point. In your case, since you do not have a global grid, you can fix the problem by setting gsnAddCyclic to False. You could instead change your code to assign the attributes lon2d and lat2d to the psfc variable rather than assigning to the resource values sfXArray and sfYArray. That seems to work okay with either setting of gsnAddCyclic. Be sure to turn off tfDoNDCOverlay whichever way you do it. This does result in a better fit of the landmask to the map outlines. Maybe the pole_lat and stand_lon values are not exactly correct.  -dave On Nov 30, 2012, at 8:12 AM, David Jones wrote: Dear Dave, > > >thank you very much for your response. I really appreciate your help. > > > >Thanks for the file you sent. I made some slight changes (i.e. exact corners, etc.) (attached) and was able to get it to work reasonably well, but it's not really completely representing the borders in places. For example, if I plot the variable LANDMASK instead of PSFC, there are gaps near the edges of, for example, the Arabian Peninsula. (Ignore the strange looking landmask in the polar regions, Greenland etc., as that is an artifact of the reanalysis data i'm using). I think that the method you've provided could well be the best that it can work using the non-wrf routines in NCL. > > > >I've uploaded the file I'm using on the ftp (it's called wrfout_d01_2009-12-01_00:00:00) so you can see this yourself if you like. I've reduced it to just one time step to make the file smaller. > > > >When I try to use the 2D arrays XLAT and XLONG to define where the map borders should be drawn I get an error message like "warning:ScalarFieldSetValues: 2d coordinate array sfXArray has an incorrect dimension size: defaulting sfXArray". Note that those 2D arrays give the geographical lat and lon for each grid point in the domain. My hope was that there would be some way to do this using the WRF routines, as WRF output is not fully NetCDF compliant - which could cause problems for some of the routines maybe. I thought that maybe this could be done through attributes like LatLonOverlay, sfXArray and sfYArray. > > > >FYI, ref_lat and ref_lon define the geographical latitude/longitude of a specific grid point in the domain. By default it is the middle of the domain, however in my case I have shifted it northwards to the grid point defined by (ref_x,ref_y) as listed in the comments at the bottom of my original script (re-attached). > > >Many thanks for your assistance, > > >Dave. > > > > >________________________________ > From: David Brown <dbrown_at_ucar.edu> >To: David Jones <jonesd647_at_yahoo.ca> >Cc: ncl-talk Talk <ncl-talk_at_ucar.edu>; wrfhelp_at_ucar.edu >Sent: Friday, November 30, 2012 5:56:38 AM >Subject: Re: Trying to get correct map over rotated lat-lon data > > >Hi David, > > >I think that rotated lat/lon projections in WRF is a confusing subject, especially in the context of NCL, because I suspect that the NCL WRF routines may not always map the attributes defined in the WRF file correctly into the resource settings required for NCL (in the rotated lat/lon case only). I will add the caveat that I do not fully understand the WRF attributes defined in your file either. What I do know is based on looking through this web page:  http://www.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm#_Description_of_GEOGRID.TBL, which may or may not be the most current description of the WRF projections and corresponding namelist parameters.  > > >Since I do not have your data file I cannot run your script for myself, but I have attempted to approximate the map projection of your data based on the WRF parameters you list in the script, the documentation listed above, and the appearance of the surface pressure field plot you attached. I am attaching my test script along with the output. I think it is fairly close to the projection of your data, but of course it would be better if I had the locations of the lower left and upper right corners of your data.  If you use gsn_csm_contour_map instead of wrf_map_overlays you should be able to plot the correct map projection for your data. One thing the documentation states is that CEN_LON and CEN_LAT do not apply for rotated lat/lon projections, but I have to wonder why they are set at all in that case. I am also not sure what the role of ret_lat and ref_lon are. I do note that ref_lat - pole_lat is almost equal to CEN_LAT and - STAND_LON is not too far
 from CEN_LON, but I don't know if that has any significance.  > > >It would be helpful if you would upload your data file. It uses a projection and a domain we have not often encountered. If you need instructions on how to do that, see http://www.ncl.ucar.edu/report_bug.shtml. Thanks, > -dave > > ><test.png> > > > > > >On Nov 28, 2012, at 11:50 AM, David Jones wrote: > > >>Dear NCL talk help, >> >> >>I would like to kindly ask for some advice regarding the plotting of rotated lat-lon output. I am currently trying to plot the output from a rotated lat-lon grid using NCL... >> >> >>When I make the contour plot of my variable (psfc), all of the contoured data appears where I expect and want it to be, the land-sea boundaries can even be detected as the value of the variable changes. So everything seems fine in that respect. >> >> >>However, when I try to plot the national/geographical borders using a map overlay, the borders are all in the wrong places. Please see the attached image, in which it is clearly visible. >> >> >>The variable I am plotting (PSFC) is made up of (time, west_east, south_north), with west_east and south_north corresponding to a point on the computational grid (i.e. they're defined by a lat and lon from the computational grid). >> >> >>I also have 2D arrays (XLAT and XLONG) giving the geographical latitude and longitude for every point on the computational grid. I tried to use these to plot the geographical boundaries in the correct location, but sadly to no avail. >> >> >>I have looked through the NCL talk archives and on the NCL website, but unfortunately nothing I've tried has worked so far. >> >> >>I'll include my script below and also attach a copy of the output from an ncl_filedump on my data. You'll notice in the script below that lots of the map options are commented out; that's because several tries showed them to make no difference whatsoever to the plot I received, whether they were included or commented out (I did try plotting with them uncommented though, in several different combinations). The only one of those options that did make a difference was setting mpres_at_mpCenterLatF = a_at_CEN_LAT. But that only changed it from one incorrect plot to another incorrect one, albeit the new looks slightly better (attached). >> >> >>If you had any guidance as to where I might be going wrong I would be most appreciative. >> >> >>Many thanks for any advice, >>Dave. >> >> >>----------------------SCRIPT----------------------------- >> >> >>; Basic script for plotting model output >> >> >>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/wrf/WRFUserARW.ncl" >>load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" >> >> >>begin >> >> >>a = addfile("wrfout_d01_2009-12-01_00:00:00.nc","r") >> >> >>wks = gsn_open_wks("x11","PSFC") >> >> >>;--------------------------------------------------------------- >>; Set basic resources >>res = True >>res_at_MainTitle = "PSFC" ; Give a main title >>pltres = True ; Plotting resources >>mpres = True ; Map resources >> >> >>;--------------------------------------------------------------- >>; Get variables >>;---1. WRF------------------ >> >> >>psfc = a->PSFC(0,:,:) >>psfc_size = dimsizes(psfc(:,:)) >>print("dimsizes(psfc) = " + psfc_size) >> >> >>;---2. Lat/Lon arrays------- >> >> >>lat2d = a->XLAT(0,:,:) >>printVarSummary(lat2d) >>lon2d = a->XLONG(0,:,:) >>printVarSummary(lon2d) >> >> >>;-------------------------------------------------------------- >>; Plotting options >>opts = res ; Add basic resources >>opts_at_FieldTitle = "PSFC" >>opts_at_cnFillOn = True ; Shaded plot >>opts_at_gsnSpreadColorEnd = -8 ; >> >> >>contour_psfc = wrf_contour(a,wks,psfc,opts) >> >> >>; Set some map options >> >> >>; mpres_at_tfDoNDCOverlay= True >> >> >>; pltres_at_sfXArray = lon2d >>; pltres_at_sfYArray = lat2d >>; pltres_at_LatLonOverlay = True >> >> >>; mpres_at_mpCenterLonF = a_at_CEN_LON ;= -16.84654 >>mpres_at_mpCenterLatF = a_at_CEN_LAT ;= 53.4331 >> >> >>; mpres_at_mpProjection = "CylindricalEquidistant" >>; mpres_at_mpLimitMode = "Corners" >>; mpres_at_mpLeftCornerLatF = lat2d(0,0) ; bottom left >>; mpres_at_mpLeftCornerLonF = lon2d(0,0) ; bottom left >>; mpres_at_mpRightCornerLatF= lat2d(268,568) ; top right >>; mpres_at_mpRightCornerLonF= lon2d(268,568) ; top right >> >> >>mpres_at_mpGeophysicalLineColor = "Black" >>mpres_at_mpGridLineColor = "Black" >>mpres_at_mpLimbLineColor = "Black" >>mpres_at_mpNationalLineColor = "Black" >>mpres_at_mpUSStateLineColor = "Black" >>mpres_at_mpPerimLineColor = "Black" >>mpres_at_mpGeophysicalLineThicknessF = 1.0 >>mpres_at_mpGridLineThicknessF = 1.0 >>mpres_at_mpLimbLineThicknessF = 1.0 >>mpres_at_mpNationalLineThicknessF = 1.0 >>mpres_at_mpUSStateLineThicknessF = 1.0 >> >> >>; Namelist map projection information: >>; map_proj : "lat-lon" (cylindrical equidistant) >>; ref_lat : 88.5 (geographical lat of grid point (ref_x,ref_y) >>; ref_lon : -5.0 (geographical lon of grid point (ref_x,ref_y) >>; ref_x : 285 (total x grid points = 570) >>; ref_y : 265 (total y grid points = 270) >>; pole_lat : 35.0 (lat of geog N hem pole on comp grid) >>; pole_lon : 180.0(lon of geog N hem pole on comp grid) >>; stand_lon: 17.5 (rotation about geographic poles) >> >> >>; MAKE PLOTS >> >> >>plot = wrf_map_overlays(a,wks,(/contour_psfc/),pltres,mpres) >> >> >>delete(res) >> >>end >> >><ncl_filedump.txt><psfc_rotlatlon.ncl><PSFC.000001.png>_______________________________________________ >>ncl-talk mailing list >>List instructions, subscriber options, unsubscribe: >>http://mailman.ucar.edu/mailman/listinfo/ncl-talk >> > > ><test_rlatlon.ncl><psfc_rotlatlon.ncl>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk