using the data output by wrf ，an exception occurs when plotting the rain distribution #167

Blissful-Jasper · 2022-03-29T01:21:41Z

Using ERA5 data as the initial field, after wrf simulation, the output file has an error in plotting the precipitation distribution, and a white line appears in the figure. But there is no problem using NCL
result and code show as below：

from wrf import getvar,get_cartopy
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter 
from netCDF4 import Dataset as nc
import cmaps

ncfile=nc(r'D:/wrf/out_24/wrfout_d01_2004-06-24_18_00_00')

slp = getvar(ncfile, "slp")
ter = getvar(ncfile, "HGT")

rain_exp=getvar(ncfile, "RAINNC")/25.4
rain_con=getvar(ncfile, "RAINC")/25.4
rain_tot=np.array(rain_exp+rain_con)*6

lons=np.array(rain_exp.XLONG)
lats=np.array(rain_exp.XLAT)

#######plot####################
cart_proj = get_cartopy(slp)
scale = '50m' 
fig = plt.figure(figsize=(12,9))
ax = plt.axes(projection=cart_proj)
land = cfeature.NaturalEarthFeature('physical',
                                        'land',
                                        scale,
                                        edgecolor='g', 
                                        facecolor='none'
                                      )
ax.add_feature(land)  # set land 
ax.coastlines(scale)
plot=ax.contourf(lons,lats,rain_tot,cmap=cmaps.CBR_coldhot,levels=np.arange(0,5,0.25),
             transform=ccrs.PlateCarree())
cb=plt.colorbar(plot, ax=ax, orientation="horizontal", pad=.05)
plt.title("Total rain")
plt.show()
####

The text was updated successfully, but these errors were encountered:

michaelavs · 2022-03-29T16:16:35Z

Hi @xpji,
I haven't been able to replicate this issue locally, but I am also using a different dataset that I have with your provided code. Is there a way I could access your dataset so I can better evaluate what the issue may be? If you would like to move this conversation to a slightly less public forum, I would recommend requesting access to the wrf-python-talk google group (I am an admin there as well) where we can then share the file via email or google drive link.

Overall, NCL does a lot of things somewhat "behind the scenes", so I assume that the issue may be something to do with how you are defining the bounds of this data projection in Python (based on levels in contourf vs the colorbar levels not being evenly spaced). If you can't share the dataset with me, we can definitely try some debugging/trouble shooting here, but the scope may be limited because of that.

Blissful-Jasper · 2022-03-30T01:35:08Z

hello，@michaelavs
Thank you for your reply, I started to find a way to upload the data to Google Cloud Drive yesterday, the link is here. If you can download this data, please check it for me, thank you.
https://drive.google.com/file/d/1dVULjkagi9VUVbFyGlXe9-QRy_tvu7mS/view?usp=sharing
At the same time，the levels in the code are set according to the example given by NCL. In addition, I found that white lines will appear when the levels are specified. If the levels are not specified, the following situation will occur.

ncfile=nc(r'D:/wrf/out_24/wrfout_d01_2004-06-24_18_00_00')

slp = getvar(ncfile, "slp")
ter = getvar(ncfile, "HGT")

rain_exp=getvar(ncfile, "RAINNC")
rain_con=getvar(ncfile, "RAINC")
rain_tot=np.array(rain_exp+rain_con)

lons=np.array(rain_exp.XLONG)
lats=np.array(rain_exp.XLAT)

#######plot####################
cart_proj = get_cartopy(slp)
scale = '50m' 
fig = plt.figure(figsize=(12,9))
ax = plt.axes(projection=cart_proj)
land = cfeature.NaturalEarthFeature('physical',
                                        'land',
                                        scale,
                                        edgecolor='g', 
                                        facecolor='none'
                                      )
ax.add_feature(land)  # set land 
ax.coastlines(scale)
plot=ax.contourf(lons,lats,rain_tot,cmap=cmaps.CBR_coldhot,
             transform=cart_proj)
cb=plt.colorbar(plot, ax=ax, orientation="horizontal", pad=.05)
plt.title("Total rain")
plt.show()

Looking forward to your reply！！

michaelavs · 2022-03-31T14:26:43Z

Hi @xpji,
Would you be able to provide the NCL code you were using that worked with this dataset? I think it would be useful for me to see the functions that were used so I can better understand what should be happening when using Python instead.

Blissful-Jasper · 2022-03-31T15:00:29Z

hi@michaelavs ，of course，the code of NCL show as below:

;*************************************************
; WRF_pcp_1.ncl
;
; Concepts illustrated:
;   - Plotting WRF data
;   - Overlaying WRF precipitation on terrain map using wrf_xxx functions
;   - Changing the size of a PNG image
;   - Creating two contour plots with two sets of filled contours
;   - Creating a color map using RGB triplets
;   - Creating a color map using RGBA quadruplets
;   - Explicitly setting contour levels
;   - Adding a title to a labelbar
;   - Drawing fully transparent filled contours
;   - Creating horizontal and vertical labelbars
;   - Adding a vertical title to a labelbar
;
; These files are loaded by default in NCL V6.2.0 and newer
; 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"

;----------------------------------------------------------------------
; This function returns an earth color palette
;----------------------------------------------------------------------
function earth_pal()
local cmap
begin
  cmap = (/ (/ 25,  25, 255/), (/ 20, 170,  42/), (/ 28, 175,  35/), \
            (/ 36, 179,  29/), (/ 44, 184,  22/), (/ 51, 189,  16/), \
            (/ 59, 194,   9/), (/ 67, 198,   2/), (/ 70, 200,   0/), \
            (/ 71, 200,   1/), (/ 72, 199,   1/), (/ 73, 199,   2/), \
            (/ 75, 198,   2/), (/ 76, 198,   3/), (/ 77, 197,   3/), \
            (/ 78, 197,   4/), (/ 79, 197,   4/), (/ 80, 196,   5/), \
            (/ 82, 196,   5/), (/ 84, 195,   6/), (/ 85, 194,   7/), \
            (/ 86, 194,   7/), (/ 87, 194,   8/), (/ 89, 193,   8/), \
            (/ 90, 193,   9/), (/ 91, 192,   9/), (/ 92, 192,  10/), \
            (/ 93, 191,  10/), (/ 94, 191,  11/), (/ 95, 191,  11/), \
            (/ 97, 190,  12/), (/ 98, 190,  12/), (/ 99, 189,  13/), \
            (/100, 189,  13/), (/101, 188,  14/), (/102, 188,  14/), \
            (/104, 188,  15/), (/105, 187,  15/), (/106, 187,  16/), \
            (/107, 186,  16/), (/108, 186,  17/), (/109, 185,  17/), \
            (/111, 185,  18/), (/112, 185,  18/), (/113, 184,  19/), \
            (/114, 184,  19/), (/115, 183,  20/), (/116, 183,  20/), \
            (/118, 182,  21/), (/119, 182,  22/), (/120, 182,  22/), \
            (/121, 181,  23/), (/122, 181,  23/), (/123, 180,  24/), \
            (/124, 180,  24/), (/126, 180,  25/), (/127, 179,  25/), \
            (/128, 179,  26/), (/129, 178,  26/), (/130, 178,  27/), \
            (/131, 177,  27/), (/133, 177,  28/), (/134, 177,  28/), \
            (/135, 176,  29/), (/136, 176,  29/), (/137, 175,  30/), \
            (/138, 175,  30/), (/140, 174,  31/), (/141, 174,  31/), \
            (/142, 174,  32/), (/143, 173,  32/), (/144, 173,  33/), \
            (/145, 172,  33/), (/146, 172,  34/), (/148, 171,  34/), \
            (/149, 171,  35/), (/150, 171,  35/), (/151, 170,  36/), \
            (/152, 170,  36/), (/153, 169,  37/), (/155, 169,  37/), \
            (/156, 168,  38/), (/157, 168,  38/), (/158, 168,  39/), \
            (/159, 167,  39/), (/160, 167,  40/), (/162, 166,  40/), \
            (/163, 166,  41/), (/164, 165,  42/), (/165, 165,  42/), \
            (/165, 165,  43/), (/165, 165,  44/), (/166, 166,  45/), \
            (/166, 166,  46/), (/166, 166,  47/), (/166, 166,  48/), \
            (/166, 166,  49/), (/167, 167,  50/), (/167, 167,  51/), \
            (/167, 167,  52/), (/167, 167,  53/), (/168, 168,  54/), \
            (/168, 168,  55/), (/168, 168,  56/), (/168, 168,  57/), \
            (/169, 169,  58/), (/169, 169,  59/), (/169, 169,  60/), \
            (/169, 169,  61/), (/169, 169,  62/), (/170, 170,  63/), \
            (/170, 170,  64/), (/170, 170,  65/), (/170, 170,  66/), \
            (/171, 171,  67/), (/171, 171,  68/), (/171, 171,  69/), \
            (/171, 171,  70/), (/171, 171,  71/), (/172, 172,  72/), \
            (/172, 172,  73/), (/172, 172,  74/), (/172, 172,  75/), \
            (/172, 172,  76/), (/173, 173,  77/), (/173, 173,  78/), \
            (/173, 173,  79/), (/173, 173,  80/), (/174, 174,  81/), \
            (/174, 174,  82/), (/174, 174,  83/), (/174, 174,  84/), \
            (/175, 175,  85/), (/175, 175,  86/), (/175, 175,  87/), \
            (/175, 175,  88/), (/175, 175,  89/), (/176, 176,  90/), \
            (/176, 176,  91/), (/176, 176,  92/), (/176, 176,  93/), \
            (/177, 177,  94/), (/177, 177,  95/), (/177, 177,  96/), \
            (/177, 177,  97/), (/177, 177,  98/), (/178, 178,  99/), \
            (/178, 178, 100/), (/178, 178, 101/), (/178, 178, 102/), \
            (/178, 178, 103/), (/179, 179, 104/), (/179, 179, 105/), \
            (/179, 179, 106/), (/179, 179, 107/), (/180, 180, 108/), \
            (/180, 180, 109/), (/180, 180, 110/), (/180, 180, 111/), \
            (/181, 181, 112/), (/181, 181, 113/), (/181, 181, 114/), \
            (/181, 181, 115/), (/181, 181, 116/), (/182, 182, 117/), \
            (/182, 182, 118/), (/182, 182, 119/), (/182, 182, 120/), \
            (/183, 183, 121/), (/183, 183, 122/), (/183, 183, 123/), \
            (/183, 183, 124/), (/183, 183, 125/), (/184, 184, 126/), \
            (/184, 184, 127/), (/184, 184, 128/), (/184, 184, 129/), \
            (/184, 184, 130/), (/185, 185, 131/), (/185, 185, 132/), \
            (/185, 185, 133/), (/185, 185, 134/), (/185, 185, 135/), \
            (/186, 186, 135/), (/186, 186, 136/), (/186, 186, 137/), \
            (/186, 186, 138/), (/187, 187, 139/), (/187, 187, 140/), \
            (/187, 187, 141/), (/187, 187, 142/), (/187, 187, 143/), \
            (/188, 188, 144/), (/188, 188, 145/), (/188, 188, 146/), \
            (/188, 188, 147/), (/188, 188, 148/), (/189, 189, 149/), \
            (/189, 189, 150/), (/189, 189, 151/), (/189, 189, 152/), \
            (/190, 190, 153/), (/190, 190, 154/), (/190, 190, 155/), \
            (/190, 190, 156/), (/190, 190, 157/), (/191, 191, 158/), \
            (/191, 191, 159/), (/191, 191, 160/), (/191, 191, 161/), \
            (/191, 191, 162/), (/192, 192, 162/), (/192, 192, 163/), \
            (/192, 192, 164/), (/192, 192, 165/), (/193, 193, 166/), \
            (/193, 193, 167/), (/193, 193, 168/), (/193, 193, 169/), \
            (/193, 193, 170/), (/194, 194, 171/), (/194, 194, 172/), \
            (/194, 194, 173/), (/194, 194, 174/), (/194, 194, 175/), \
            (/195, 195, 176/), (/195, 195, 177/), (/195, 195, 178/), \
            (/195, 195, 179/), (/196, 196, 180/), (/196, 196, 181/), \
            (/196, 196, 182/), (/196, 196, 183/), (/196, 196, 184/), \
            (/197, 197, 185/), (/197, 197, 186/), (/197, 197, 187/), \
            (/197, 197, 188/), (/197, 197, 189/), (/198, 198, 189/), \
            (/198, 198, 190/), (/198, 198, 191/), (/198, 198, 192/), \
            (/199, 199, 193/), (/199, 199, 194/), (/199, 199, 195/), \
            (/199, 199, 196/), (/199, 199, 197/), (/200, 200, 198/), \
            (/200, 200, 199/), (/200, 200, 200/) /)/255.
  return(cmap)
end

;----------------------------------------------------------------------
; This function returns a precipitation color palette
;----------------------------------------------------------------------
function precip_pal()
local cmap, ncolors
begin
  cmap = (/ (/242, 242, 242/), (/154, 192, 205/), (/178, 223, 238/), \
            (/191, 239, 255/), (/  0, 235, 235/), (/  0, 163, 247/), \
            (/  0, 255,   0/), (/  0, 199,   0/), (/  0, 143,   0/), \
            (/  0,  63,   0/), (/255, 255,   0/), (/255, 143,   0/), \
            (/255,   0,   0/), (/215,   0,   0/), (/191,   0,   0/), \
            (/255,   0, 255/), (/155,  87, 203/), (/ 92,  52, 176/) /)

;---Create RGBA array from RGB array so we can set a couple colors to transparent 
   ncolors          = dimsizes(cmap(:,0))
   cmap_rgba        = new((/ncolors,4/),float)
   cmap_rgba(:,0:2) = cmap/255.
   cmap_rgba(:,3)   = 1.0     ; set all colors to fully opaque
   cmap_rgba(0:1,3) = 0.0     ; set first two colors to fully transparent
   return(cmap_rgba)
end


;----------------------------------------------------------------------
; Main code
;----------------------------------------------------------------------
begin
;
; Open file and get variables.
;
  a = addfile("/out_24/wrfout_d01_2004-06-24_18_00_00","r")

;
; Terrain
;
  ter = wrf_user_getvar(a,"HGT",0)
  ter@description = "Terrain Height"
  ter@units       = "m"
;
; Get non-convective, convective
; Calculate total precipitation
;
  scale = 6.       ; Artificial scale factor to see all colors
  rain_exp = wrf_user_getvar(a,"RAINNC",-1) /25.4
  rain_con = wrf_user_getvar(a,"RAINC", -1) /25.4
  rain_tot = (rain_exp + rain_con)*scale
  rain_tot@description = "Total Precipitation (inches)"

; Define contour levels here so we can determine up front
; what the merged color maps should be.

  precip_levels = (/ .01, 0.25, 0.50, 0.75, \
                    1.00, 1.25, 1.50, 1.75, \
                    2.00, 2.25, 2.50, 2.75, \
                    3.00, 3.25, 3.50, 3.75, \
                    4.00, 4.25, 4.50, 4.75, \
                    5.00 /)

  ter_levels = (/ .01, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, \
                   15, 16 /) * 100.

  wtype          = "png"
; wtype@wkWidth  = 1500     ; Increase size for a slightly
;  wtype@wkHeight = 1500     ; better looking PNG.
  wks = gsn_open_wks(wtype,"WRFpcp")    ; Send graphics to PNG file

;
; Set up resource list that will be shared between the
; two wrf_contour calls.
;
  res                      = True
  res@gsnDraw              = False
  res@gsnFrame             = False
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnFillOn             = True
  res@cnLinesOn            = False
;
; Generate plot of terrain for background.
;
; First set up resource list specific to terrain plot.
;
  opts_ter                          = res

  opts_ter@cnLevels                 = ter_levels
;
; In order to use cnFillPalette with wrf_contour, you must
; gsnSpreadColors to False and cnSpanFillPalette to True.
;
  opts_ter@gsnSpreadColors          = False
  opts_ter@cnSpanFillPalette        = True
  opts_ter@cnFillPalette            = earth_pal()

; Resources to control labelbar title, size, and location.
  opts_ter@lbTitleString            = "Terrain"
  opts_ter@pmLabelBarWidthF         = 0.8
  opts_ter@pmLabelBarHeightF        = 0.35
  opts_ter@pmLabelBarOrthogonalPosF = -0.18

; Create terrain contours (no drawing done yet).
  contour_ter = wrf_contour(a,wks,ter,opts_ter)

; Plotting options for precipitation
  opts_r                      = res

  opts_r@cnLevels             = precip_levels

  opts_r@gsnSpreadColors          = False
  opts_r@cnSpanFillPalette        = True
  opts_r@cnFillPalette            = precip_pal()

  opts_r@cnSmoothingOn        = True
  opts_r@cnSmoothingDistanceF = .005

; Resources to control precipitation labelbar, which will be 
; vertical. 
  opts_r@lbTitleString            = "Total Precipitation"
  opts_r@lbTitleDirection         = "Down"
  opts_r@lbTitleJust              = "CenterRight"
  opts_r@lbTitlePosition          = "Right"
  opts_r@lbTitleOffsetF           = 0.07

  opts_r@lbOrientation            = "Vertical"
  opts_r@pmLabelBarSide           = "Right"
  opts_r@pmLabelBarHeightF        = 0.77
  opts_r@pmLabelBarWidthF         = 0.11
  opts_r@pmLabelBarOrthogonalPosF = 0.03
  opts_r@lbBoxMinorExtentF        = 0.4

; Total Precipitation (color fill)
  contour_tot = wrf_contour(a, wks, rain_tot(0,:,:), opts_r)

;
; Use the special wrf_map_overlays function to figure out
; the correct map projection and do the overlay.
;
; Set up two resource lists for wrf_map_overlays. 
;
  pltres                            = True
  mpres                             = True
  mpres@mpGeophysicalLineColor      = "Black"
  mpres@mpNationalLineColor         = "Black"
  mpres@mpUSStateLineColor          = "Black"
  mpres@mpGridLineColor             = "Black"
  mpres@mpLimbLineColor             = "Black"
  mpres@mpPerimLineColor            = "Black"
;  mpres@mpGeophysicalLineThicknessF = .5
;  mpres@mpUSStateLineThicknessF     = .5
;  mpres@mpDataBaseVersion           = "HighRes"
;  mpres@mpDataResolution            = "FinestResolution"


  plot = wrf_map_overlays(a,wks,(/contour_ter,contour_tot/),pltres,mpres)

end

Blissful-Jasper · 2022-04-05T03:13:18Z

Hello,@michaelavs， have you found the cause of the above problem? I use other wrf output files to find that I still get the same problem. I don't know if it is the problem of the drawing code or the file. Look forward to your reply.

michaelavs · 2022-04-05T17:24:14Z

Hi @xpji,
I've been looking through everything the past few days and have a working version of your code, but I think there are a couple of factors leading to the output you are getting. Before diving more into the working code, I'll point out what I think the issues are and why:

The wrfout file you are using is a netcdf file with HDF5 formatting. I believe that the HDF5 formatting is causing an issue with wrf-python, which appears to only accept pure netcdf files and doesn't have support for HDF5 formatting (based on the documentation page where only NetCDF is mentioned). I think this isn't a limitation in NCL (again based on documentation) and that is why the plot is fine when you run the data using NCL (which touches on the "behind the scenes" topic I mentioned in my first response). If you wanted to use the python package H5PY, that may help to resolve the issue but I was not able to fully explore this option because of what I found with the next point
Once removing all level specifications from contourf and the colorbar calls, the lines through the plot disappeared. By doing this, I let python make the contour levels based on the rain data and was able to get a better idea of what it was looking at in terms of the range of the data points. There does seem to be a weird occurrence when you exceed a certain number of levels and/or spacing of values which will then reintroduce the lines through the plot, but the location and size of them varies depending on the level values/spacing. In the code below, the values for "levels" in contourf and "ticks" in colorbar are the highest I could get before the lines would start to appear again.

Overall, I think the biggest culprit is that the data file is HDF5 formatted and that is causing some confusion for wrf-python and matplotlib overall. Another aspect to back this up is that I ran your code using a two different wrfout files, one being a different time step of the datafile used in the NCL-WRF example you were following, and did not have a problem with plotting the data, however, both of those wrfout files were pure netCDF files (i.e. the formatting wasn't HDF5, but netcdf). If this was an issue with wrf-python overall, I would not expect to have been able to successfully plot the data using your code, and should have seen the same issues we are seeing when using your dataset.

Here is the working code I was able to create:

from wrf import getvar,get_cartopy, latlon_coords
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

from netCDF4 import Dataset as nc

import cmaps


ncfile = nc("/Users/misi1684/Downloads/wrfout_d01_2004-06-24_18_00_00",
            decode_times=True)

scale = 6

rain_exp = getvar(ncfile, "RAINNC")/25.4
rain_con = getvar(ncfile, "RAINC")/25.4
rain_tot = (rain_exp+rain_con)*scale

slp = getvar(ncfile, "slp")
ter = getvar(ncfile, "HGT")

lats, lons = latlon_coords(slp)

################### plot ###################

cart_proj = get_cartopy(slp)

scale = '50m'
fig = plt.figure(figsize=(12,6))
ax = plt.axes(projection=cart_proj)
ax.coastlines(scale)

land = cfeature.NaturalEarthFeature('physical',
                                    'land',
                                    scale,
                                    edgecolor='g',
                                    facecolor='none'
                                    )
ax.add_feature(land)  # set land

plot = ax.contourf(lons,
                 lats,
                 rain_tot,
                 cmap=cmaps.CBR_coldhot,
                 levels=np.linspace(0.0, 120, 16),
                 transform=ccrs.PlateCarree(),
                 colorbar=False)

plt.colorbar(plot, ax=ax,
                    orientation="horizontal",
                    ticks=np.linspace(0.0, 120, 16),
                    drawedges=True,
                    extendrect=True,
                    pad=0.08,
                    shrink=0.75,
                    aspect=30)

plt.title("Total rain")
plt.show()

Output projection:

You can definitely play around with the levels/ticks parameters, but from what I have seen, this is the most resolved version of the dataset that does not cause the lines through the projection.

Blissful-Jasper · 2022-04-06T01:28:37Z

hi，@michaelavs
nice to receive your reply message.

The format of this file is netcdf when I choose to output, so I don't know why it is mixed with hdf5 instead of a “pure” nc format data.
I tried to use h5py package for data reading and drawing, however, the same problem occurred.
My strangest thing is why NCL draws the result so well.And why is the variable for rain in the NCL code divided by 25.4, if python does not handle the rain data like this,its optimum range and maximum value may different, there will still be a big problem,it's very confusing to me. NCL results and manipulation of data works like magic

In any case, thank you very much for communicating with you and wishing you a happy life.

michaelavs closed this as completed Apr 11, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

using the data output by wrf ，an exception occurs when plotting the rain distribution #167

using the data output by wrf ，an exception occurs when plotting the rain distribution #167

Blissful-Jasper commented Mar 29, 2022

michaelavs commented Mar 29, 2022

Blissful-Jasper commented Mar 30, 2022

michaelavs commented Mar 31, 2022

Blissful-Jasper commented Mar 31, 2022

Blissful-Jasper commented Apr 5, 2022

michaelavs commented Apr 5, 2022

Blissful-Jasper commented Apr 6, 2022

using the data output by wrf ，an exception occurs when plotting the rain distribution #167

using the data output by wrf ，an exception occurs when plotting the rain distribution #167

Comments

Blissful-Jasper commented Mar 29, 2022

michaelavs commented Mar 29, 2022

Blissful-Jasper commented Mar 30, 2022

michaelavs commented Mar 31, 2022

Blissful-Jasper commented Mar 31, 2022

Blissful-Jasper commented Apr 5, 2022

michaelavs commented Apr 5, 2022

Blissful-Jasper commented Apr 6, 2022