VERDI shapefile export - grid cells shifted to west

[lizadams]

Describe the bug Exporting LWMASK from 4km domain centered over Long Island shows that the grid cells are shifted to the West and South

To Reproduce Steps to reproduce the behavior:

1. Download 4km file containing LWMASK variable GRIDCRO2D_04OTC2_160930.nc

2. Load file into verdi and create tile plot of the LWMASK variable using command

   verdi.sh -f $cwd/GRIDCRO2D_04OTC2_160930.nc -s "LWMASK[1]" -g tile

3. Download the high resolution background map provided in USA_States.zip

4. Click on GIS Layers> Add Map Layers > Other

5. Navigate to the folder where the USA_States.shp file is located and double click on the USA_States.shp and click open.

6. Zoom in on Block Island, Rhode Island (Island located North East of Long Island) by specifing the columns and rows in the tile plot. For the VERDI tile plot use the following settings for the Rows and Columns: Click on Controls > Set Row and Column > Rows: 140 .. 144, Columns 110 .. 115

7. Find the Lat/Lon Value of the upper left grid cell where the LWMASK value is 1

   Click on Controls > Show Lat/Lon Hover the mouse over the upper left corner of the grid cell that is colored Red or has a value of 1 and record the lat/lon value displayed in the lower right corner of the Tile Plot 41.1763N, 71.59667W

8. Export the Shapefile by clicking on

   File > Export as Image/GIS > Change Files of Type to Shapefile, and set Filename to record what is being saved LWMASK_4km_GRIDCRO2D_04OTC2_160930_zoomed_in_Block_Island I uploaded a tar.gz file that contains the shapefiles: LWMASK_4km_GRIDCRO2D_04OTC2_160930_zoomed_in_Block_Island.tar.gz

9. Load the high resolution US_STATES.shp and exported 4km_LWMASK file into QGIS and check the lat/lon value of the grid cell that has LWMASK=1 near Channel Island. VERDI reports the lat/lon value as 41.1763N, 71.59667W

   

10. The value in QGIS for the upper left corner of the grid cell near Block Island is 41.1863, 71.6661

    

Expected behavior The shapefile export should have a grid cell with LWMASK=1 over Block Island when loaded to QGIS.

Version Information :

• OS: IOS -VERDI 2.0 2021-02-17

Link to google drive with relevant files. https://drive.google.com/drive/folders/1gdI8-_KN1hi-Y9WAhPnUgaBqCRkig2r3?usp=sharing

Header information in the GRIDCRO2D_04OTC2_160930.nc file can be examined using the command: ncdump -h GRIDCRO2D_04OTC2_160930.nc > GRIDCRO2D_04OTC2_160930.txt :GDTYP = 2 ; :P_ALP = 33. ; :P_BET = 45. ; :P_GAM = -97. ; :XCENT = -97. ; :YCENT = 40. ; :XORIG = 1644000. ; :YORIG = -144000. ; :XCELL = 4000. ; :YCELL = 4000. ; :VGTYP = 7 ; :VGTOP = 5000.f ; :VGLVLS = 1.f, 0.9975f ; :GDNAM = "04OTC2_CROSS " ; :UPNAM = "OUTGM3IO " ;

This file contains the LAT/LON Values, So, I should be able to extract only the grid cell that is located over Block Island, and see if the LAT/LON values in the dataset match what VERDI is showing in the Tile Plot, and also what is being exported. LAT = 37.17028, 37.16284, 37.15538, 37.14792, 37.14042, 37.13292, 37.12539, 37.11785, 37.11029, 37.10272, 37.09512, 37.08752, 37.07989, 37.07225, 37.06459, 37.05691, 37.0492, 37.04149, 37.03376, 37.02602, 37.01826, 37.01048, 37.00267, 36.99486, 36.98703, 36.97918, 36.97131, 36.96342, 36.95552, 36.9476, 36.93967, 36.93172, 36.92375, 36.91576, 36.90775, 36.89973, 36.8917, 36.88364, 36.87557, 36.86749, 36.85937, 36.85126, 36.84311, 36.83496, 36.82678, 36.81859, 36.81037, 36.80215, 36.79391,

LON = -78.19421, -78.14981, -78.10541, -78.06104, -78.01666, -77.97229, -77.92795, -77.88361, -77.83926, -77.79492, -77.75061, -77.7063, -77.66199, -77.61771, -77.57343, -77.52914, -77.48486, -77.44061, -77.39639, -77.35214, -77.30792, -77.2637, -77.21948, -77.17529, -77.1311, -77.08694, -77.04279, -76.9986, -76.95447, -76.91034,

To verify the LAT/LON values in the GRIDCRO2D_04OTC2_160930.nc match the values obtained by using Show Lat/Lon Values in VERDI, I loaded the LAT and LON values, and then used the mouse to show the values for each in a csv format.

Using the Probe, and putting it in the center of each grid, the values seem to match what is reported in the VERDI Tile Probe that was exported. Lower Grid Cell Center:

Upper Grid Cell Center:

[lizadams]

LWMASK_4km_GRIDCRO2D_04OTC2_160930_zoomed_in_Block_Island.prj contains the following information:

PROJCS["", GEOGCS["Sphere", DATUM["WGS_1984", SPHEROID["WGS 84", 6378137.0, 298.257223563], TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]], PROJECTION["Lambert_Conformal_Conic_2SP"], PARAMETER["central_meridian", -97.0], PARAMETER["latitude_of_origin", 40.0], PARAMETER["standard_parallel_1", 45.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["standard_parallel_2", 33.0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH]]

[lizadams]

Pehaps this reference will help. https://fabienmaussion.info/2018/01/06/wrf-projection/ "- if WRF were to use a datum shift to interpret the land-surface data into its spherical coordinates (as suggested by the authors of the papers), then we should also do it when mapping its results

• since WRF doesn't, we also shouldn't"

This reference refers to https://ral.ucar.edu/sites/default/files/public/WRFHydro_GIS_Preprocessor_v5.pdf

Failing to apply an appropriate transformation between spherical and spheroidal datums, such as simply treating geodetic latitudes as geocentric, can result in significant (up to 20km North-South) locational shift in the spatial features contained in the data. This shift is greatest at 45° latitude (David et al 2009). No predefined transformation exists in ArcGIS between a sphere and spheroid, thus a custom transformation must be created. As suggested by Cedric David, a ‘geocentric translation’ may be defined in ArcGIS, leaving all parameters ‘0’ when translating between sphere and spheroid

https://journals.ametsoc.org/view/journals/mwre/141/6/mwr-d-12-00351.1.xml Overlapping Interests: The Impact of Geographic Coordinate Assumptions on Limited-Area Atmospheric Model Simulations

[lizadams]

Confirmed that the VERDI plot matches the NCL plot of the GRIDCRO2D_04OTC2_160930.nc file, and that the LWMASK covers Block Island. So the LWMASK correctly identifies land in the area of Block island in both the VERDI Tile Plot, and the NCL Plot.

The script that I used to creat the ncl plot is available in the following google drive: https://drive.google.com/drive/folders/10VtTFZhHyuRHGNq_xDTYsXMgRHSLn0N9?usp=sharing To run the script, I first load the NCAR Graphics module

module load ncarg Then source the environment variable file source ncl_env.csh Then run ncl to create the plot ncl plot_LWMASK.ncl The following plot is generated. NCL_Plot_LWMASK4km_GRIDCRO2D_04OTC2_160930.png

[lizadams]

Extracted the LANDMASK variable from the wrf file, so that I would have a smaller file to work with.

module load nco ncks -v LANDMASK wrfout_d01_2016-09-29_00:00:00 wrfout_d01_2016-09-29_00:00:00.LANDMASK Loaded the wrfout file into QGIS, goal is to use QGIS to generate a shapefile. The option to Load Vector Layer with the netCDF plugin didn't work...

I found these instructions on how to create a shapefile from a wrf output file. https://gis.stackexchange.com/questions/88940/conversion-of-latitude-longitude-to-different-map-projection/123237#123237 Another option is to use WRF-python https://wrf-python.readthedocs.io/en/latest/plot.html

[yadongxuEPA]

Retested with VERDI_2.0_linux64_20210427.tar.gz, found that the projection string in the exported shapefiles was changed to the followings: PROJCS["unknown", GEOGCS["GCS_unknown", DATUM["D_unknown", SPHEROID["unknown", 6370000.0, 0.0]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Longitude", EAST], AXIS["Latitude", NORTH]], PROJECTION["Lambert_Conformal_Conic"], PARAMETER["central_meridian", -97.0], PARAMETER["latitude_of_origin", 40.0], PARAMETER["standard_parallel_1", 45.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["standard_parallel_2", 33.0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH]]

Loaded the exported shapefiles and "USA_States.shp" into QGIS to verify (as shown in the following map).
the exported grid cell was not shifted to west anymore and the coordinates for the upper left corner ([41.17766,-71.59912]) were very close to the ones shown in VERDI.

[yadongxuEPA]

Retested with VERDI_2.1_linux64_20210526.tar.gz, found that the projection string in the exported shapefiles was changed to the followings: PROJCS["unknown", GEOGCS["GCS_unknown", DATUM["D_unknown", SPHEROID["unknown", 6370000.0, 0.0]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Longitude", EAST], AXIS["Latitude", NORTH]], PROJECTION["Lambert_Conformal_Conic"], PARAMETER["central_meridian", -97.0], PARAMETER["latitude_of_origin", 40.0], PARAMETER["standard_parallel_1", 45.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["standard_parallel_2", 33.0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH]] Loaded the exported shapefiles and "USA_States.shp" into QGIS to verify (as shown in the following map). the exported grid cell was not shifted to west anymore. This issue has resolved.

[yadongxuEPA]

Retested with VERDI_2.1_linux64_20210626.tar.gz, found that the projection string in the exported shapefiles was changed to the followings: PROJCS["unknown", GEOGCS["GCS_unknown", DATUM["D_unknown", SPHEROID["unknown", 6370000.0, 0.0]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH]], PROJECTION["Lambert_Conformal_Conic_2SP"], PARAMETER["central_meridian", -97.0], PARAMETER["latitude_of_origin", 40.0], PARAMETER["standard_parallel_1", 45.0], PARAMETER["false_easting", 0.0], PARAMETER["false_northing", 0.0], PARAMETER["scale_factor", 1.0], PARAMETER["standard_parallel_2", 33.0], UNIT["m", 1.0], AXIS["x", EAST], AXIS["y", NORTH]] Loaded the exported shapefiles and "USA_States.shp" into QGIS to verify (as shown in the following map). Although the projection string has changed slightly from the previous one, the location of the exported grid cell is accurate.

[yadongxuEPA]

Retested with VERDI_2.1_linux64_20210906.tar.gz, confirmed that the upgraded openJDK has not caused any problem on this issue.