VERDI Tileplot Display of LANDMASK variable from geogrid polar stereographic dataset is incorrect

[lizadams]

Describe the bug Meghan Mallard created a polar stereographic geogrid file that is 4km resolution. She also provided a lambert conformal geogrid file that is of the same 4km domain. VERDI is able to accurately plot the LANDMASK variable for the lambert conformal 4km geogrid file, but not for the polar stereographic 4km geogrid file.

To Reproduce Steps to reproduce the behavior:

1. Loaded the files into verdi and created tile plots of the LANDMASK variable.
2. ./verdi.sh -f $cwd/../geogrid_polar_stereographic/geo_em.polar_stereo_d03.nc -s "LANDMASK[1]" -g tile -f $cwd/../geogrid_lambert_conformal/geo_em.lambert_conformal_d03.nc -s "LANDMASK[2]" -g tile
3. Used the Externalize view into floating window, and then repositioned the two plots so they could be viewed side by side.

Expected behavior Expect that these plots both have the landmask over the same area, and follow the coastlines of the detailed map. There are several differences 1) The Lambert Conformal dataset is displayed using the North America Background Map by default. 2) The Polar Stereographic dataset is displayed using the World Map by default. 3) When I switch to using the same North Americal Background Map for both plots, it looks like the scale factor is different, that the lambert conformal data is "zoomed in" compared to the polar stereographic data that is displayed.

4) The lambert conformal dataset display doesn't cover as far North and West as the polar stereographic dataset. (You can see Mexico in the polar stereographic dataset, but not in the lambert conformal dataset).

Using verdi build 3-27-2021

Datasets are available from this Google Drive:

https://drive.google.com/drive/u/0/folders/1-Rd497PLYuG6NJuCBnQtA1A7RkvExp

Note, I renamed the original files so we can keep track of what projection they are in when visualizing them in VERDI.

There are 3 different geogrid files for the polar stereographic projection 36km geo_em.polar_stereo_d01.nc 12km geo_em.polar_stereo_d02.nc 4km geo_em.polar_stereo_d01.nc

There are 3 different geogrid files for the lambert conformal projection 36km: geo_em.lambert_conformal_d01.nc 12km: geo_em.lambert_conformal_d02.nc 4km: geo_em.lambert_conformal_d03.nc

The lambert conformal projection files are visualized correctly in VERDI The polar stereographic projection files are not visualized correctly in VERDI.

I am going to try to do a shapefile export of the data, just to see what it looks like when I bring the data into QGIS.

Liz

[lizadams]

When I try to export the two datasets as shapefiles and then visualize them in VERDI.

I get what looks like an incomplete shapefile for the polar stereographic dataset. 100 Apr 2 17:28 landmask_geo_em.polar_stereo_d03.shx -rw-rw-r-- 1 lizadams rc_cep-emc_psx 100 Apr 2 17:28 landmask_geo_em.polar_stereo_d03.shp -rw-rw-r-- 1 lizadams rc_cep-emc_psx 543 Apr 2 17:28 landmask_geo_em.polar_stereo_d03.prj -rw-rw-r-- 1 lizadams rc_cep-emc_psx 65 Apr 2 17:28 landmask_geo_em.polar_stereo_d03.dbf -rw-rw-r-- 1 lizadams rc_cep-emc_psx 632 Apr 2 17:30 landmask_geo_em.lambert_conformal_d03.prj -rw-rw-r-- 1 lizadams rc_cep-emc_psx 37362700 Apr 2 17:30 landmask_geo_em.lambert_conformal_d03.shp -rw-rw-r-- 1 lizadams rc_cep-emc_psx 6868190 Apr 2 17:30 landmask_geo_em.lambert_conformal_d03.dbf -rw-rw-r-- 1 lizadams rc_cep-emc_psx 3296713 Apr 2 17:30 landmask_geo_em.lambert_conformal_d03.fix -rw-rw-r-- 1 lizadams rc_cep-emc_psx 2197900 Apr 2 17:30 landmask_geo_em.lambert_conformal_d03.shx

I am able to correctly visualize the lambert_conformal shapefile after I edit the datum from wgs1984 to a sphere. I can't visualize the polar stereographic shapefile, all of the LANDMASK variable values are zero.

[yadongxuEPA]

Retested with VERDI_2.1_linux64_20210626.tar.gz, found that the issue is not resolved. The polar stereographic projection files are still not visualized correctly in VERDI (shown as below) and the exported shapefiles were incomplete and could not be visualized in QGIS.

[yadongxuEPA]

Retested with VERDI_2.1_linux64_20210707.tar.gz, found that the issue is not resolved. The polar stereographic projection files are still not visualized correctly in VERDI (shown as below).

[yadongxuEPA]

Retested with VERDI_2.1_linux64_20210809.tar.gz, found that the issue is not resolved. The polar stereographic projection files are still not visualized correctly in VERDI (shown as below).

[lizadams]

I confirmed that NCARG also plots the polar sterographic plot incorrectly.

module load ncarg

Generated the plots using the following command: ncl wrf_metgrid2.ncl

cat wrf_metgrid2.ncl

;   Example script to plot skintemp from a single metgrid file
;   November 2008

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
  a = addfile("./geo_em.polar_stereo_d01.nc","r")  ; Open a file

; We generate plots, but what kind do we prefer?
  type = "x11"
; type = "pdf"
; type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"plt_metgrid_2")

  res = True                           ; Set up some basic plot resources
  res@MainTitle = "METGRID FILES geo_em.polar_stereo_d01.nc"
  res@Footer = False

  pltres = True
  mpres = True

    skint = wrf_user_getvar(a,"LANDMASK",0)    ; Get skintemp from file

    opts = res
    opts@cnFillOn            = True
    contour = wrf_contour(a,wks,skint,opts)
    plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)

end

[lizadams]

I tried another NCL script to see if that would work any better, and it has the same issues. ;***** ; WRF: color over map with lat/lon labels ;**** 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/wrf/WRF_contributed.ncl" begin ;**** ; open file and read in data ;**** f = addfile ("geo_em.lambert_conformal_d01.nc", "r") ;**** ; Read character variable Times; Convert to string for plots ; Read vertical coordinate for plot labels ;**** times = chartostring(f->Times) ; built-in function ;znu = f->ZNU(0,:) ; (Time, bottom_top) ;**** ; Read perturbation geopotential at all times and levels ;**** x = f->GREENFRAC ; (Time, bottom_top, south_north, west_east) ;**** ; create plots ;**** wks = gsn_open_wks("ps" ,"geo_em.lambert_conformal_d01_GREENFRAC") ; ps,pdf,x11,ncgm,eps gsn_define_colormap(wks,"BlAqGrYeOrReVi200") ; select color map

res = True ; plot mods desired res@gsnMaximize = True ; uncomment to maximize size res@gsnSpreadColors = True ; use full range of colormap res@cnFillOn = True ; color plot desired res@cnLinesOn = False ; turn off contour lines res@cnLineLabelsOn = False ; turn off contour labels

;**** ; Use WRF_contributed procedure to set map resources ;**** ; WRF_map_c (f, res, 0) res = wrf_map_resources(f,res) ;**** ;**** ; set True for native projection (faster) ;**** res@tfDoNDCOverlay = True

;**** ; Turn on lat / lon labeling ;**** res@pmTickMarkDisplayMode = "Always" ; turn on tickmarks

;**** ; Loop over all times and levels ( uncomment ) ; Demo: one arbitrarily closen time and level ;**** dimx = dimsizes(x) ; dimensions of x ntim = dimx(0) ; number of time steps klev = dimx(1) ; number of "bottom_top" levels

nt = ntim/2 ; arbitrary time kl = 6 ; " level ;;do nt=0,ntim-1 ; uncomment for loop ;; do ll=0,klev-1 res@tiMainString = times(nt) res@gsnLeftString = x@description+" z=" plot = gsn_csm_contour_map(wks,x(nt,kl,:,:),res) ;; end do ;;end do end

geo_em.lambert_conformal_d01_GREENFRAC.pdf geo_em.polar_stereo_d01_GREENFRAC.pdf

[lizadams]

I was able to get wrf-python to plot the geo_em.polar_stereo_d01.nc file correctly. I had to reinstall anaconda3, create an environment called geo_env, and then use conda-forge I followed the instructions on how to create the environment and use it on the geopandas website. https://geopandas.org/getting_started/install.html

conda install -c conda-forge cartopy Eventually, I was able to follow the instructions on the wrf-python site to generate the following plot using jupyter notebook.

https://wrf-python.readthedocs.io/en/latest/plot.html#matplotlib-with-cartopy

from netCDF4 import Dataset
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as crs
from cartopy.feature import NaturalEarthFeature

from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim,
                 cartopy_ylim, latlon_coords)

# Open the NetCDF file
#ncfile = Dataset("/Users/lizadams/Downloads/wrf_python_tutorial/wrf_tutorial_data/wrfout_d01_2005-08-28_00_00_00")
ncfile = Dataset("/Users/lizadams/Downloads/geogrid_polar_stereographic/geo_em.polar_stereo_d01.nc")
# Get the sea level pressure
slp = getvar(ncfile, "LANDMASK")

# Smooth the sea level pressure since it tends to be noisy near the
# mountains
smooth_slp = smooth2d(slp, 3, cenweight=4)

# Get the latitude and longitude points
lats, lons = latlon_coords(slp)

# Get the cartopy mapping object
cart_proj = get_cartopy(slp)

# Create a figure
fig = plt.figure(figsize=(12,6))
# Set the GeoAxes to the projection used by WRF
ax = plt.axes(projection=cart_proj)

# Download and add the states and coastlines
states = NaturalEarthFeature(category="cultural", scale="50m",
                             facecolor="none",
                             name="admin_1_states_provinces_shp")
ax.add_feature(states, linewidth=.5, edgecolor="black")
ax.coastlines('50m', linewidth=0.8)

# Make the contour outlines and filled contours for the smoothed sea level
# pressure.
plt.contour(to_np(lons), to_np(lats), to_np(smooth_slp), 10, colors="black",
            transform=crs.PlateCarree())
plt.contourf(to_np(lons), to_np(lats), to_np(smooth_slp), 10,
             transform=crs.PlateCarree(),
             cmap=get_cmap("jet"))

# Add a color bar
plt.colorbar(ax=ax, shrink=.98)

# Set the map bounds
ax.set_xlim(cartopy_xlim(smooth_slp))
ax.set_ylim(cartopy_ylim(smooth_slp))

# Add the gridlines
ax.gridlines(color="black", linestyle="dotted")

plt.title("LANDMASK geo_em.polar_stereo_d01.nc")

plt.show()

[lizadams]

I used the same jupyter notebook to plot the geo_em file that was created using the lambert_conformal projection.

geo_em.lambert_conformal_d01.nc

[lizadams]

So, VERDI is not able to read and use the projectin information correctly for the polar stereographic geo_em.nc file. I am including a link that may help Tony diagnose the issue. https://fabienmaussion.info/2018/01/06/wrf-projection/

Specifically, this information: (note this example is using lambert conformal conic projection, and the settings necessary to specify the polar stereographic projection may vary slightly. In any case, we need to examine how VERDI reads the dataset projection information from the WRF input file and sets the projection information. This is a different routine than how VERDI reads information from an MCIP or CMAQ file.

The projection parameters are stored as attributes in the NetCDF file:

import pyproj wrf_proj = pyproj.Proj(proj='lcc', # projection type: Lambert Conformal Conic lat_1=ds.TRUELAT1, lat_2=ds.TRUELAT2, # Cone intersects with the sphere lat_0=ds.MOAD_CEN_LAT, lon_0=ds.STAND_LON, # Center point a=6370000, b=6370000) # This is it! The Earth is a perfect sphere.

We need to use something similar in VERDI, by reading in the projection information from the WRF file.

[lizadams]

We have different types of input files, and VERDI needs to use some method to recognize them.

The geo*.nc files are part of WRF, and use the wrf options to specify the map projection, see this PDF for more information: http://140.112.69.65/research/coawst/COAWST_TUTORIAL/training_2019/monday/werner_wps.pdf

The GRIDCRO2D files are I/O API based, and in VERDI use the M3IO Convention, with the GDTYP variable to specify the map projection.

The MPAS files are also treated as a separate convention.

I used wrf-python to correctly plot both the polar stereographic and the lambert conformal versions of the WRF geogrid file.

geogrid_polar_stereographic/geo_em.polar_stereo_d01.nc and geogrid_lambert_conformal/geo_em.lambert_conformal_d01.nc

NCVIEW can plot the data, but when you use the OPT button to add the coastline, the data is not correctly georeferenced for either geo_em* file. I has also tried using NCVIEW with other WRFOUT files, and they also fail to correctly georeference once you add the coastline. So, I don't think we can gleam any information from the source code for NCVIEW. On dogwood: module load ncview ncview ./geogrid_lambert_conformal/geo_em.lambert_conformal_d01.nc

Click on 2D vars, select LANDMASK variable Click on OPT and then select USA or Coastline

I think the section of the VERDI code that needs to be fixed is here: https://github.com/CEMPD/VERDI/blob/9108c508d0b75c43ee813daba805a568f27f00a9/verdi_data_loaders/src/ucar/nc2/dataset/conv/WRFConvention.java#L83

Starting at line 271.

Comparing what is done for M3IO Convention
https://github.com/CEMPD/VERDI/blob/9108c508d0b75c43ee813daba805a568f27f00a9/verdi_data_loaders/src/ucar/nc2/dataset/conv/M3IOConvention.java#L404

M3IO Convention starting on line 390 for Polar Stereographic:

private CoordinateTransform makePolarStereographicProjection(NetcdfDataset ds) {
double central_meridian = findAttributeDouble(ds, "P_GAM");
double xcent = findAttributeDouble(ds, "XCENT");
double ycent = findAttributeDouble(ds, "YCENT");
// double par1 = findAttributeDouble(ds, "P_ALP");
double par2 = findAttributeDouble(ds, "P_BET");

// Stereographic sg_tmp = new Stereographic(ycent, xcent, 1.0, 0, 0, earthRadius);
// sg_tmp.setCentralMeridian(central_meridian);
// double calc_scale = sg_tmp.getScale();

double scale2 = (1.0 + Math.sin(Math.toRadians(par2))) / 2.0;
// System.out.println("in Stereographic: calculated scale2 " + scale2);

Stereographic sg = new Stereographic(ycent, xcent, scale2, 0, 0, earthRadius);
sg.setCentralMeridian(central_meridian);
return new ProjectionCT("PolarStereographic", "FGDC", sg);
}

WRF Convention (case 2 is polar stereographic)

else {
double lat1 = findAttributeDouble(ds, "TRUELAT1");
double lat2 = findAttributeDouble(ds, "TRUELAT2");
double centralLat = findAttributeDouble(ds, "CEN_LAT"); // center of grid
double centralLon = findAttributeDouble(ds, "CEN_LON"); // center of grid
double standardLon = findAttributeDouble(ds, "STAND_LON"); // true longitude
double standardLat = findAttributeDouble(ds, "MOAD_CEN_LAT");
ProjectionImpl proj = null;
switch (projType) {
case 0: // for diagnostic runs with no georeferencing
proj = new FlatEarth();
projCT = new ProjectionCT("flat_earth", "FGDC", proj);
// Logger.debug(" using LC "+proj.paramsToString());
break;
case 1:

proj = new LambertConformal(standardLat, standardLon, lat1, lat2, 0.0, 0.0, 6370.);
projCT = new ProjectionCT("Lambert", "FGDC", proj);

break;
case 2:
// Thanks to Heiko Klein for figuring out WRF Stereographic
double lon0 = (Double.isNaN(standardLon)) ? centralLon : standardLon;
double lat0 = (Double.isNaN(centralLat)) ? lat2 : centralLat; // ?? 7/20/2010
double scaleFactor = (1 + Math.abs( Math.sin(Math.toRadians(lat1)))) / 2.; // R Schmunk 9/10/07
// proj = new Stereographic(lat2, lon0, scaleFactor);
proj = new Stereographic(lat0, lon0, scaleFactor);

Can we replace this with something that is similar to what we use for the M3IO convention proj = new Stereographic(lat0, lon0, scaleFactor); ^^^ replace with: proj = new Stereographic(lat0, lon0, scaleFactor , 0, 0, earthRadius);

proj.setCentralMeridian(central_meridian);

^^ Do we need to set the Central Meridian for the Polar Stereographic projection? I have commented it out above as I am not sure. return new ProjectionCT("PolarStereographic", "FGDC", proj);

^^^^^

projCT = new ProjectionCT("Stereographic", "FGDC", proj); break;

[lizadams]

I confirmed that IDV also has the same error when plotting the polar stereographic geo_em* files.

[lizadams]

I am trying to use gdalinfo

gdalinfo geo_em.polar_stereo_d01.nc > gdalinfo.geo_em.polar_stereo_d01 gdalinfo geo_em.lambert_conformal_d01.nc > gdalinfo.geo_em.lambert_conformal_d01

Then take a difference of the gdalinfo output for the two files diff ../geogrid_polar_stereographic/gdalinfo.geo_em.polar_stereo_d01 ./gdalinfo.geo_em.lambert_conformal_d01

2c2
< Files: geo_em.polar_stereo_d01.nc
---
> Files: geo_em.lambert_conformal_d01.nc
6,9c6,9
<   NC_GLOBAL#CEN_LAT=36.549969
<   NC_GLOBAL#CEN_LON=-91.370155
<   NC_GLOBAL#corner_lats={8.8275099,47.072079,38.778862,6.0575075,8.7806873,46.919132,38.626015,6.002058,8.7160788,47.16972,38.854023,5.9557781,8.669425,47.016289,38.700821,5.900516}
<   NC_GLOBAL#corner_lons={-119.73759,-154.36772,-33.251747,-68.448944,-119.85035,-154.51091,-33.155357,-68.346695,-119.6904,-154.59279,-33.055305,-68.504547,-119.80296,-154.73546,-32.959427,-68.402412}
---
>   NC_GLOBAL#CEN_LAT=36.774815
>   NC_GLOBAL#CEN_LON=-91.51355
>   NC_GLOBAL#corner_lats={8.3103256,53.973969,48.97438,5.8370514,8.2693863,53.890453,48.87307,5.7864838,8.1739197,54.106762,49.098553,5.707222,8.1330795,54.023098,48.997063,5.6567917}
>   NC_GLOBAL#corner_lons={-123.43092,-147.89371,-34.908508,-63.314636,-123.5688,-148.11934,-34.719513,-63.184143,-123.38971,-148.03575,-34.754089,-63.365295,-123.52737,-148.2617,-34.564941,-63.234955}
28c28
<   NC_GLOBAL#MAP_PROJ=2
---
>   NC_GLOBAL#MAP_PROJ=1
30c30
<   NC_GLOBAL#MOAD_CEN_LAT=36.549969
---
>   NC_GLOBAL#MOAD_CEN_LAT=36.774815
54c54
<   SUBDATASET_1_NAME=NETCDF:"geo_em.polar_stereo_d01.nc":Times
---
>   SUBDATASET_1_NAME=NETCDF:"geo_em.lambert_conformal_d01.nc":Times

Is the information from gdalinfo for these files consistent with what you are using in VERDI?

[lizadams]

I've scoured AWS public data, but can't find WRF output that was created using the Polar Stereographic projection.

This thread mentions the issue that we are seeing. https://www.ncl.ucar.edu/Support/talk_archives/2010/0381.html

and they mention a cause of the error being https://www.ncl.ucar.edu/Support/talk_archives/2010/0377.html

The issue is old, and the links don't work.

OK - I just found an wrfout file that is polar stereographic from a WRF-CHEM tutorial. https://ruc.noaa.gov/wrf/wrf-chem/wrf_tutorial_emissions_v35/exercise_4.html

Found a wrfoutput file that is polar stereographic

VERDI gives an error message when trying to plot it, and it seems to be shifted in a similar way to the geo_em file. wrfout_d01_polar_stereographic.nc.tar.gz

./verdi.sh -f /Users/lizadams/downloads/wrfout_d01_polar_stereographic.nc 2021.08.21 21:14:58.646 [main] ERROR anl.verdi.loaders.NetcdfBoxer - Error while creating CRS for Stereographic Error in "GEOGCS": Can't parse "PRIMEM["Greenwich",0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS (...) 0], PARAMETER["false_northing", 0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH]]" because "PRIMEM" is unrecognized.

The best variable to plot to compare to the western US coastline was the LU_INDEX.

[lizadams]

Also tried using several other visualization packages to view both the geo_em and the WRF output data with the Polar Stereographic Projection.

Salem https://salem.readthedocs.io/en/v0.3.4/wrf.html

Both Salem and WRF-Python are able to plot the geo_em files correctly for both the lambert conformal and the polar stereographic projection.

IDV tutorial for WRF https://www.eoas.ubc.ca/courses/atsc507/A507Resources/WRF_running/wrf-tutorial3_idv_display.pdf

NCL https://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/NCL/Examples/BASIC/wrf_single_field.php

WRF-Python Tutorial https://wrf-python.readthedocs.io/en/latest/tutorials/wrf_workshop_2019.html

[lizadams]

When I use the python package salem to display the WRF output file that I found that had the polar stereographic projection, it gave a warning, and didn't seem to display the LU_INDEX. It did show the map in the same general location that VERDI gave.

[lizadams]

At this time, I do not recommend providing WRF (geo_em or wrfout) files in the data/model directory with the Polar Stereographic Projection. We should make a notation in the release notes that VERDI has a bug when displaying WRF files with the Stereographic Projection.

[yadongxuEPA]

Retested this issue with VERDI_2.1_linux64_20211206.tar.gz, the issue has resolved. A tile plot for the "LANDMASK" variable in "geo_em.polar_stereo_d03.nc" was correctly displayed in VERDI main window, shown as below: Then the exported shapefiles were correctly displayed in QGIS, shown as below:

[lizadams]

I also tested the WRF file wrfout_d01_polar_stereographic.nc.tar.gz using the LU_INDEX variable, and it was correctly displayed in VERDI.