using the data output by wrf ，an exception occurs when plotting the rain distribution
Closed this issue · 7 comments
Blissful-Jasper commented
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()
####
michaelavs commented
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 commented
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 commented
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 commented
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 commented
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 commented
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 commented
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.