Upside down & unknown projection/CRS & large-value backwater of streamflow in CHRTOUT_GRID1 output

[Martin Nguyen]

Hi team,

Thank you very much for the great WRF-Hydro modelling system. I am new to this. I am trying to analyse the streamflow for my case. However, the streamflow in my CHRTOUT_GRID1 seems to have three issues on my side:

- It's upside down

- There is not much information about its projection/crs. I tried to changed it to 'Lambert Conformal Conic' which is NZGD2000/NZCS2000 EPSG:3851 on my side, but still not succeed.

- It has a lot of backwater with very large values. I found some relevant questions here WRF-Hydro negative streamflow. However, I did not understand much about the timestep in the answer.

Could you please help me with these issues? Please let me know if you need further information. I really appreciate!

[aubrey]

Hi Martin:

How are you creating your geogrid and routing stack? Are you using our tools? Are you using one of our supported projections? And when you say your CHRTOUT_GRID is upside down, what tool are you using to view this file? Does it look upside down in a simple netcdf viewer like ncview, or are you trying to import it into a GIS?

Re timesteps, you can check our user guide for guidance on how to set timesteps based on your grid spacing: https://wrf-hydro.readthedocs.io/en/latest/index.html

Hope that helps.

Aubrey

[Martin Nguyen]

Hi Aubrey,

Thank you so much for your email. I really apprecite it. I used your tool to create geogrid and routing stack as below:

- For geogrid: I used this command <./geogrid/geogrid.exe>. The namelist.wps I attached as below. I tried 'lambert' and 'lat-lon' but they both produces the same results as I described above. Also, if I used matplotlib and xarray to plot the streamflow and topography in python, they are on the same scale (meters) but the streamflow is smaller than the DEM.

- For routing stack: I used your python code from here  wrf_hydro_gis_preprocessor/wrfhydro_gis at master · NCAR/wrf_hydro_gis_preprocessor. The code is <python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_020_full_100m\geo_em.d01.nc --CSV C:\Users\mng42\wrf_wps\example_nz_020_full_100m\river_source.csv -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_020_full_100m\dem_mataura_largescale_filled.tif -R 5 -t 1500 -o C:\Users\mng42\wrf_wps\example_nz_020_full_100m\output\mataura_test_001.zip>

For the 'supported projections', I'm not so sure about this, so far as I know, I think I did not use because I don't know about that. Could you please share a bit more details about this?

For the tool I used to view it, I used QGIS to view it, I set the openstreetmap at EPSG:4326, the topography (DEM or Fullhires_dom.TOPOGRAPHY) is aligning with the openstreetmap, but the streamflow is as I described above. When I used ncview, it did not show anything, only white blank screen appeared.

For importing into a GIS, I'm not so sure about this, I think I did not try to do it. However, could you please give me a bit more details about this?

I have attached here two namelist.wps I used - one is lambert and the other is latlon. Please let me know if you would like me to provide the DEM or the results. I can send you through your email.

Kind regards

Martin Nguyen 

[Aubrey Dugger]

Hi Martin:

You can find specifications on the projections supported in the GIS pre-processor in the user guide: https://github.com/NCAR/wrf_hydro_arcgis_preprocessor/releases/download/5.2.0/WRFHydro_GIS_Preprocessor_v5_2_0.pdf

I'm not well versed in WPS, so can't comment on your setup there. But I know we support lambert conformal conic as I use that regularly. My guess is that the CHRTOUT_GRID file is not fully CF compliant and QGIS is not interpreting the y coordinates. You can probably fix this by swapping the order of the y coordinates (I think WPS/WRF/Noah-MP use south-to-north, while a lot of GIS applications use north-to-south). There is a simple NCO command to do this:

ncpdq -O -a -y filename.nc filename.nc

I can view CHRTOUT_GRID files OK in ncview.

You might try pulling time series from CHRTOUT_DOMAIN files. If you want to share your model setup I can take a look. As always, we recommend starting your testing with idealized forcings to make sure your domain and setup are working.

Thanks!

Aubrey

--

-----------------------------------------------------------

Aubrey Dugger

NCAR Research Applications Laboratory

Office: 303-497-8418, Cell: 310-663-5115

adu...@ucar.edu

[Kevin Sampson]

Martin,

The outputs from WPS for a Lambert Conformal Conic domain and a Cylindrical Equidistant (lat-lon) domain will not be identical. I can see how they might seem that way when looking at the domain using xarray, because the dimensions of the domains specified in your wps namelist files are identical. However, these will indeed be written to different coordinate systems. Review the projection information in the global attributes of the resulting geo_em.d01.nc file for those differences. You can read more about these projections here.

I suggest that you review your domain by converting the routing stack to GeoTiff using the WRF-Hydro GIS Pre-processor (Examine_Outputs_of_GIS_Preprocessor.py). This will automatically interpret the coordinate system, and output properly formatted GeoTiff files that you can view and overlay in GIS.

Alternatively, you can obtain the definition of the coordinate system of your domain by outputting a GIS projection file (.prj) using the WRF-Hydro GIS Pre-processor tool "Build_PRJ_From_Geogrid_File.py".

Keep in mind that WPS writes output files from south-to-north with the y-coordinate ascending. The WRF-Hydro GIS Pre-processing system, for historical reasons, writes the Fulldom_hires.nc file from north-to-south, consistent with GDAL-derived datasets. This could be the reason you see the streamflow outputs as being upside-down, if rendering with a library that does not properly interpret the cf-metadata.

Thanks,

Kevin

--
You received this message because you are subscribed to the Google Groups "wrf-hydro_users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wrf-hydro_use...@ucar.edu.
To view this discussion visit https://groups.google.com/a/ucar.edu/d/msgid/wrf-hydro_users/CA%2BOs1ReqaWBEAvVY%2B0R%3Ds1_V822AWVZT7dqwbfoTsi%3DaxOWHGA%40mail.gmail.com.

[Martin Nguyen]

Hi Audrey and Kevin,

Thank you very much for your all your replies. I have been trying as you suggested but I still faced those problems to read all of them in QGIS. However, when I read the TOPOGRAPHY and streamflow from CHRTOUT_GRID1 on a same plot in python using xarray, they seem like on a same projection, but it looks like the WRF-Hydro only produce the streamflow (in CHROUT_GRID1) for a part of the domain, not full domain (as attached picture). Also, the values of streamflow are from negative (backwater) to really high positive values. Could you please help me with these issues? 

I attached here my namelist.wps, my domain (Fulldom_hires.nc, etc.), Forcing data (from WPS and regridded already), and other necessary data. Please let me know if you cannot access it. I also attached here the plot I used in python. You can see the streamflow is not predicted thoroughly for the full domain. 

Building stack code: <python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_020_full_100m\geo_em.d01.nc --CSV C:\Users\mng42\wrf_wps\example_nz_020_full_100m\river_source.csv -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_020_full_100m\dem_mataura_largescale_filled.tif -R 5 -t 1500 -o C:\Users\mng42\wrf_wps\example_nz_020_full_100m\output\mataura_test_001.zip>

Kind regards

Martin Nguyen

[Aubrey Dugger]

Hi Martin:

Have you tested your model with idealized forcings (FORC_TYP=4)? Do you still see the strange streamflow values? 

Thanks!

Aubrey

[Martin Nguyen]

Hi Aubrey,

Thank you so much for your reply. I chose FORC_TYP = 1. I have tried FORC_TYP = 4 and the problem is the same. However, thanks to this, I have noticed that our forcing data is actually the WRF output, so we have the forcing input variable names as below:

! Forcing input variable names

forcing_name_T = "T2"

forcing_name_Q = "Q2"

forcing_name_U = "U10"

forcing_name_V = "V10"

forcing_name_P = "PSFC"

forcing_name_LW = "GLW"

forcing_name_SW = "SWDOWN"

forcing_name_PR = "RAINNC"

In that the RAINNC has unit of mm NOT mm/s or kg/m^2 like RAINRATE as in this instruction. I think this might be the problem that causes the strangely large and small values in streamflow. I have tried to change into FORC_TYP = 3 (which is WRF or 7 which is for WRF w/ spec. precip.), but our Forcing data is stored in Windows folder, so the filename does not allow colon ":" as instructed in that instruction, so the WRF-Hydro could not find/read the forcing data. As we aim to store these large forcing data in windows, do you have any suggestions for this?

Thank you very much

Kind regards

Martin Nguyen

[Aubrey Dugger]

Hi Martin:

Yes, if you are using wrfouts as forcing files you should be using FORC_TYP=3, which will handle the variable names and unit conversions needed. This assumes your WRF-Hydro geogrid file is exactly the same dimensions and resolution as your WRF setup that generated the wrfouts.

I don't have a solution for the Windows issue except that you could probably edit the code that constructs the filename string for the reads:

https://github.com/NCAR/wrf_hydro_nwm_public/blob/7ec01b41ec7a8a8067a5cc82f8e1d395ab201f84/src/Routing/module_lsm_forcing.F90#L985

and then edit the filenames to remove the colons.

Alternatively you could convert the wrfout files to LDASIN format (no colons) and manually do the needed forcing variable name changes and unit conversions.

Let us know how it goes.

Thanks!

Aubrey

[Martin Nguyen]

Hi Aubrey,

Thank you so much for your reply. I have tried the way to recompile after edit the code as you pointed out. Specifically, I changed (as below) from lines 986-1007:

########################

     if(FORC_TYP.eq.3) then

!!Create forcing data filename...

        inflnm = trim(indir)//"/"//&

             "wrfout_d0"//hgrid//"_"//&

             olddate(1:4)//"-"//olddate(6:7)//"-"//olddate(9:10)//&

             "_"//olddate(12:13)//"_00_00"

        inquire (file=trim(inflnm), exist=fexist)

        if ( .not. fexist ) then

           print*, "no forcing data found", inflnm

           call hydro_stop("In read_hydro_forcing_seq() - no forcing data found")

        endif

        do i_forcing = 1, int(24*3600/dt)

           wrf_dt = i_forcing*dt

           call geth_newdate(out_date,olddate,nint(wrf_dt))

           inflnm2 = trim(indir)//"/"//&

             "wrfout_d0"//hgrid//"_"//&

             out_date(1:4)//"-"//out_date(6:7)//"-"//out_date(9:10)//&

             "_"//out_date(12:13)//"_00_00"

           inquire (file=trim(inflnm2), exist=fexist)

           if (fexist ) goto 991

#####################

After that I recompile and tried again. I changed FORC_TYP=3 to use the WRF outputs. The filename of forcing data is like this, for example, "wrfout_d03_2020-02-01_00_00_00". However, I can still not produce the simulations. I attached the run.log for reference. It looks like wrf-hydro cannot find the files in the forcing folder. Do you have any idea about this? Please let me know if you need further information.

[Martin Nguyen]

Hi Aubrey,

Just in case, the previous message did not reach you. I organised my updated message here again:

Thank you very much for your instruction.

I chose the way to rewrite the code and recompile, because our project prefers using the results directly produced from the WRF. I changed the code as you mentioned here. Specifically, at line 991 and line 1005, I changed ":00:00" into "_00_00" and then recompiled. I also changed FORC_TYP=3 to use the WRF outputs.  The filename of forcing data now looks like this "wrfout_d03_2020-02-01_00_00_00". However, I can still not produce the simulations (but if I changed back to FORC_TYPE=1, and changed the name of forcing data accordingly, it can produce simulations but with strange values of streamflow as I described before). I attached the run.log reference here. The final error message is as below, and it looks like wrf-hydro can still not find the files in the forcing folder to me. Could you please help me with this? Please let me know if you need further information.

### Final error #######

[C002FV:16183] 1 more process has sent help message help-mpi-api.txt / mpi-abort
[C002FV:16183] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages

Kind regards

Martin Nguyen

[Kevin Sampson]

  Hi Martin,

I wanted to provide a few notes after examining your domain:

1. I see that you selected the cylindrical equidistant projection (MAP_PROJ=6). Given the small area domain and the mid-latitude location, I would select Lambert Conformal Conic (MAP_PROJ=1). Lambert is appropriate for mid-latitude domains, and has been tested extensively with WRF-Hydro tools.
2. Your domain does not include the full watershed for the river system that you are trying to model. I would not expect streamflow or other processes to be well represented without the full contributing watershed.
3. I also noted that the CHRTOUT_DOMAIN1 files contain correctly-placed river points, but only for a subset of your gridded domain. I have never before seen this behavior, and I wonder if the choice of projection is the reason for this.
4. For some reason, the projection information was not accurately written to your output files (missing the latitude of origin), causing the grids to be located on the equator when exported to GeoTiff. I will look into this further, but I suspect if you change to Lambert Conformal Conic, this issue will go away.

Thanks,

Kevin

To view this discussion visit https://groups.google.com/a/ucar.edu/d/msgid/wrf-hydro_users/12c47c79-0928-49c0-acbb-f626b4cba4efn%40ucar.edu.

[Aubrey Dugger]

Hi Martin:

Can you share a few of your sample wrfout files?

Thanks!

Aubrey

[Martin Nguyen]

Hi Kevin,

Thank you so much for your email. I really appreciate it. I am trying your suggestion and will keep you updated soon.

Hi Aubrey,

Thank you so much for your email. I have placed my "FORCING_fromWRFoutputs" folder here (please let me know if you cannot access). It includes two subfolders - before and after regridding. For the regridding, I followed the instruction here (WRF* not for window machine). Please let me know if you need further information.

Kind regards

Martin Nguyen

[Martin Nguyen]

Hi Kevin,

I am trying to redo the work as you suggested. I increase the domain so it can include the source of upstream (watershed). The DEMs that I have here are only in "EPSG: 4326" or "EPSG: 2193". I tried to reproject it to the Lambert Conformal Conic (using the crs from wrf-hydro example Croton). The xmin, ymin, xmax, ymax = -20368581.4167000018060207,-1962432.7912000000942498 , -19697023.7393000014126301,-1366172.2441000000108033. 

I am building up the routing stack using the python work here. However, I have some errors while building it up and I can't find the solution. Could you please help me with this? I attached here the namelist.wps, error, domain shapefile, and the code. The DEM is quite large, so I cannot provide here. Please let me know if you need further information. I will try to find a way to send to you.

Building routing stack code: <python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_030_full_100m\geo_em.d01.nc --CSV C:\Users\mng42\wrf_wps\example_nz_030_full_100m\mataura_frxst_pts_csv.csv -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_030_full_100m\mataura_002.tif -R 10 -t 2000 -o C:\Users\mng42\wrf_wps\example_nz_030_full_100m\output\mataura_test_001.zip>

Thank you very much

Kind regards

Martin Nguyen

[Aubrey Dugger]

Hi Martin:

One note: you will want to customize the lambert conformal conic projection parameters for your specific location. I'm not sure where your basin of interest is (New Zealand?), but the Croton test case projection specifications are from a contiguous U.S. model domain. They may not work well for your region of interest, and definitely not for a whole different hemisphere. I would look for an example closer to your area.

See: https://proj.org/en/stable/operations/projections/lcc.html

Aubrey

[Martin Nguyen]

Hi Aubrey,

Thank you very much for your note. I have tried to convert it again as you suggested, unfortunately I still got the same error (as in the log file) when I tried to build up routing stack. Could you please help with this? I really apologise if my case is a bit confusing. Please let me know if you need further information. My steps to build up routing stack are as below:

STEP 1. I converted my DEM in New Zealand from EPSG:2193 to Lambert Conformal Conic using projection as the format I found on the website you sent: "+proj=lcc +lat_1=-46.5 +lat_2=-45.5 +lat_0=-46.0 +lon_0=169.5 +R=6370000 +x_0=0 +y_0=0 +units=m +no_defs". My xmin, xmax, ymin, ymax are as below:

xmin = -20368581.4167
xmax = -19697023.7393
ymin = -1962432.7912
ymax = -1366172.2442

STEP 2. After that, I did some calculation to obtain the center point of the projection-converted DEM and convert it to degrees (EPSG:4326) to input to the namelist.wps (I tried to input the Lambert Conformal Conic format but yielded errors) - I attached my namelist.wps in this email.

STEP 3. Then I created:

* <geo_em.d01.nc> using the code <./geogrid/geogrid.exe>

* <soil_properties.nc> using the code <Rscript create_soilproperties.R --geogrid="geo_em.d01.nc">

* <wrfinput_d01.nc> using the code <Rscript create_wrfinput.R --geogrid="geo_em.d01.nc">

STEP 4. Then I tried to get the domain to check using this code <python Create_Domain_Boundary_Shapefile.py -i C:\Users\mng42\wrf_wps\example_nz_030_full_100m\geo_em.d01.nc -o C:\Users\mng42\wrf_wps\example_nz_030_full_100m\output>. The xmin, xmax, ymin, ymax of the domain are as below (they are different from what I expected, very different from the ones in the step 1) - I attached the domain in this email:

xmin = -58903.6407,

xmax = 58896.35922,

ymin = -78900.5166,

ymax = 78899.48334

STEP 5. After that I tried with the "building routing stack" to see how it goes using this code <python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_030_full_100m\geo_em.d01.nc --CSV C:\Users\mng42\wrf_wps\example_nz_030_full_100m\mataura_frxst_pts_csv.csv -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_030_full_100m\mataura_002.tif -R 10 -t 500 -o C:\Users\mng42\wrf_wps\example_nz_030_full_100m\output\mataura_test_001.zip> and I still got the same error - I attached the domain in the log file in this email.

The error is RuntimeError: C:\Users\mng42\wrf_wps\example_nz_030_full_100m\output\scratchdir\streams.tif, band 1: Failed to compute statistics, no valid pixels found in sampling. I got this error at  Process: CHANNELGRID written to output netCDF. More information is in the log file I attached.

Kind regards

Martin Nguyen

[Kevin Sampson]

Martin,

Thank you for the update. 

No need to reproject your source DEM; the WRF-Hydro GIS Pre-processing tools will do this automatically, while handling the datum transformation correctly, for you. You can use your source DEM, as long as the elevation units are in meters.

I took a look at the error you are getting, and I think it has to do with your input CSV file, which I do not have and cannot examine. Can you try without the --CSV option? That way, we might be able to get past this error (for now), and see how the rest of the domain looks.

Kevin

[Martin Nguyen]

Hi Kevin and Aubrey,

Thank you so much for your instruction. I have tried without CSV and could past the error and produced the outputs. I attached here a piece of time of CHRTOUT_GRID1 for streamflow and run log (successful one) for checking. I still have some issues:

- Upside down streamflow in QGIS, I guess it was because QGIS reads the coordinates differently as you and Aubrey helped me explained before. I will try to find more solutions as my project requires it should be read by QGIS. Please let me know if you have any solutions.

- The streamflow and TOPOGRAPHY data are not aligned in QGIS (I used the one from Fulldom_hires as attached in mataura_test_001.zip)

- When I tried on python they all align well and not upside down at all. I will use python for visualisation for temporary. Please let me know if you have any other suggestions for QGIS.

- Large values and negative values for streamflow. I guess this is because I use the forcing data from WRF model output (I changed the name into *.LDASIN_DOMAIN1, FORC_TYP=1, to be able to run the wrf-hydro), the RAINNC (unit is in mm, not mm/s or kg/m^2) is the reason making the streamflow values strange like that. If using FORC_TYP = 3, I don't know why the wrf-hydro cannot work (based on the error, it looks like the wrf-hydro cannot find out my forcing data files). I attached the WRF output after regridding (please note it is under 7zip file) and the run log error (failed one) for further reference.

For the csv file for forecast points, I tried it again with changing the points into different places. It worked this time, but the frxst_pts.tif did not show the forecast point using the <Build_Routing_Stack.py>, but when I checked by using the <Examine_Outputs_of_GIS_Preprocessor.py>, the result "frxst_pts.tif" did show it. Please help me check on this one.

The code fore Build_Routing_Stack.py is:  python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_033_full_100m\geo_em.d01.nc -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_033_full_100m\mataura_4326.tif -R 4 -t 2000 -o C:\Users\mng42\wrf_wps\example_nz_033_full_100m\output\mataura_test_001.zip

Apart from those files attached, I also included other files - namelist.wps, hydro.namelist, and namelist.hrldas. Please let me know if you need further information.

I really appreciate for all your patience and big helps.

Kind regards

Martin Nguyen

[aubrey]

Hi Martin:

My guess is something is going wrong in your forcing processing. Unfortunately you cannot just rename the wrfout files to LDASIN since the variables expected from each are different. For example, wrfout precipitation is an accumulation which the wrf-hydro code adjusts internally to a rate. If you want to process from wrfout to LDASIN, you'll need to add these variable name changes and units remapping to your workflow.

I also noticed the wrfout file you sent is 399x799, but the original domain you shared is 348x898. These need to be the same dimensions (or maybe you changed the wrf-hydro setup along the way). I can test the model read of your wrfout files under FORC_TYP = 3, but I need a consistent setup to test.

Let us know how it goes.

Thanks!

Aubrey

[Martin Nguyen]

Hi Aubrey,

Thank you very much for your email. As the files are large to be attached, I have shared all the necessary files to run a simulation in vers_002_consistentsetup. I also converted the LDASIN back to the original form which is the WRF model output and provided the "before-after" regridding results. I also provided the errors (results) I got when running simulation with this version under 'Errors_of_running_simulation'. Please note that I did not include the forecast points (csv) in this simulation version, but I attached it anyway, just in case.

Please let me know if you cannot access to the folder I shared.

Building Routing stack:  <python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_040_mataura\geo_em.d01.nc -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_040_mataura\mataura_4326.tif -R 4 -t 2000 -o C:\Users\mng42\wrf_wps\example_nz_040_mataura\output\mataura_test_001.zip>

Kind regards

Martin Nguyen

[Martin Nguyen]

Hi Aubrey, Yongxin, and Kevin,

Thank you very much for your helps and suggestions. I really appreciate that. I have successfully run the simulation with WRF outputs as FORCING. The reason I could not run before was that I missed some lines of codes that I did not change the ":" to "_" in the wrfoutput filename here. It works when I changed all of them and recompiled! So now the streamflow in CHRTOUT_GRID1 is in m3s-1 and have no large or backwater values anymore.

However, I can still not overlay this streamflow on the topography (DEM) because they don't have the CRS or on Openstreetmap in QGIS. Here is how I understand how to use the wrf-hydro projection. Please help me check if I understand them correctly.

1. Since my DEM is in EPSG:4326 (lat and lon are degree) so I just need to find out the centroid of this raster for the ref_lon and ref_lat (they require the degree) for the namelist.wps. I will then keep the truelat1 and truelat2 = ref_lat, and stand_lon =  ref_lat. I choose map_proj = "lambert". Then create the geo_em.d01.nc

2. After creating the geo_dem.d01.nc and run the WRF-Hydro to get the streamflow in CHRTOUT_GRID1. The shape of streamflow at this time is really large, upside down, and has no crs and locates at some random place in QGIS.

3. I then use Build_PRJ_From_Geogrid_File.py to create the projection and write this projection into the streamflow of CHRTOUT_GRID1. The streamflow is not upside down anymore but its shape is really small and stays within the DEM.

4. After that writing that projection to the streamflow, I then convert this streamflow (the one that I just write the projection to) to the new projection (EPSG:4326) and try overlaying it on DEM/Openstreetmap using QGIS (because the DEM or the TOPOGRAPHY from Build_Routing_Stack.py or from Examine_Outputs_of_GIS_Preprocessor.py overlays each other well and DEM has crs of EPSG:4326). The problem is the same as step 3: The streamflow is not upside down anymore but its shape is really small and stays within the DEM and overlay with the one from step 3.

Could you please help me with this issue? I attached here the namelist.wps, geo_em.d01.prj (from Build_PRJ_From_Geogrid_File.py) the streamflow I extracted from CHRTOUT_GRID1, wrote projection (from Build_PRJ_From_Geogrid_File.py), and converted to 4326. For all the necessary files to run simulation and the simulation results (including the CHROUT_GRID1), I put them here - vers_003_streamflow. 

Please let me know if you need further information.

Kind regards

Martin Nguyen

[Kevin Sampson]

Martin,

I have a few questions and comments related to the issues that you are facing. 

1. Did you provide a spatial metadata file in the LAND_SPATIAL_META_FLNM parameter in your hydro namelist? This file comes out of the GIS Pre-processing tools and provides WRF-Hydro with the CF-compliant spatial metadata that is machine readable by many GIS applications and libraries. 

2. The CHRTOUT_GRID appears upside down depending on how you read the data. The routing grids in WRF-Hydro are typically written in a north-south direction due to compatibility with GIS applications. The LSM grid is written south-north. This is important information when reading and displaying the files. If you wish to overlay the files and you are reading the data without a geospatially-aware application or library, then this grid will appear upside down. You can always flip it for display.

3. It looks to me as though your streamflow data is not being placed randomly. It is likely located at the center of your projected coordinate system and assumed to have horizontal resolution in the units of the projeciton, which is likely meters. This indicates to me that the projeciton is being read or imposed properly, but the horizontal resolution of the grid is not.

4. Can you please provide your geogrid file (geo_em.d01.nc) as well as a CHRTOUT file so that I can test the coordinate system information and metadata? Ideally, if you provide the full domain directory, namelists, and all WRF-Hydro output files from the last simulation timestep, this should give me what I need to perform some tests. 

Thanks,

Kevin

To view this discussion visit https://groups.google.com/a/ucar.edu/d/msgid/wrf-hydro_users/CA%2BSiuXKTYr8ABzwLHBqUN609Ad2bPYLkQKHPzq%3Dgm4G5_t%2Ba8Q%40mail.gmail.com.

[Martin Nguyen]

Hi Kevin,

Thank you so much for your reply.

1. Yes, I provided the LAND_SPATIAL_META_FLNM in the hydro.namelist.

2 and 3.  When you said "The LSM grid is written south-north", do you mean (for example) the TOPOGRAPHY in Fulldom_hires.nc is written south-north. If so, I checked this TOPOGRAPHY and it overlays on the Openstreetmap of QGIS well. As I checked the streamflow of CHRTOUT_GRID1, it looks like the CRS is correctly set (Lambert Conformal Conic), and the files includes a GeoTransform attribute indicating the spatial resolution (50-meter grid). Thanks to your idea in 2 and 3, to ensure proper alignment with the DEM (Fuldom_hires.nc), I used the python code below to reproject the streamflow of CHRTOUT_GRID1. It does overlay the TOPOGRAPHY and the Openstreetmap in QGIS. Could you please let me know if this is the correct way?

### Python code to convert streamflow projection ################

# Load DEM

dem = rxr.open_rasterio(r"S:\FloodRiskResearch\Martin\WRF-Hydro\simulation_forKevin\DOMAIN\Fulldom_hires.nc", mask_and_scale=True)

dem_single = dem.squeeze()  # Remove extra dimensions if needed

# Load CHRTOUT streamflow

ds = xr.open_dataset(r"S:\FloodRiskResearch\Martin\WRF-Hydro\simulation_forKevin\202002010300.CHRTOUT_GRID1", mask_and_scale=True)

streamflow = ds.streamflow.squeeze()

# Flip streamflow vertically (Y is reversed in some CHRTOUT exports)

streamflow_flipped = streamflow.isel(y=slice(None, None, -1))

# Attach correct CRS and transform (from DEM) — even though CHRTOUT has it in attributes

streamflow_flipped.rio.write_crs(dem_single.rio.crs, inplace=True)

streamflow_flipped.rio.write_transform(dem_single.rio.transform(), inplace=True)

# Save as GeoTIFF

streamflow_flipped.rio.to_raster(r"S:\FloodRiskResearch\Martin\WRF-Hydro\reproject\streamflow_new_004.tif")

######################################################################################

4. As the files to provide is too large, I put them into this shared folder - vers_003_streamflow. Please let me know if you cannot access to it. It includes two subfolders and a file:

- Necessary inputs: includes namelist.wps, and results from geogrid.exe, soilproperties and wrfinput_d01 from R script, zip file from Build_Routing_Stack.py (stored in output folder), and other files

- Simulation: includes DOMAIN, FORCING, hydro.namelist, namelist.hrldas, outputs (including CHRTOUT_GRID files), and other files.

- Flipped_streamflow.tiff: from the code above (using 202002010300.CHRTOUT_GRID1)

Code for Build_Routing_Stack.py: python Build_Routing_Stack.py -i C:\Users\mng42\wrf_wps\example_nz_forKevin_001\geo_em.d01.nc --CSV C:\Users\mng42\wrf_wps\example_nz_forKevin\mataura_frxst_pts_csv.csv -b False -r True -d C:\Users\mng42\wrf_wps\example_nz_forKevin_001\mataura_4326_merged_fill.tiff -R 4 -t 2000 -o C:\Users\mng42\wrf_wps\example_nz_forKevin_001\output\mataura_test_001.zip

5. The frxst_pts in Fulldom_hires.nc does not include forecast points that I stated in the CSV file, but it does appear in the Examine_Outputs_of_GIS_Preprocessor.py. I provided the results of these python codes in the output folder of necessary_inputs (as I said above). Could you please help me check this? Those python codes I downloaded from here.

I am using WRF-Hydro version 5.4.0. Please let me know if you could not access the folder or need further information.

Kind regards

Martin Nguyen

To view this discussion visit https://groups.google.com/a/ucar.edu/d/msgid/wrf-hydro_users/CAP82fYYUmB6tbu8iC4yoB2UErqiFieWHGFybMhRnF0S-dEXUkw%40mail.gmail.com.

[Kevin Sampson]

Martin,

Thank you for the detailed explanation and for sending the requested files. 

When I stated that the LSM grid is written in a south-north orientation, I mean that the array is written to disk with the southernmost row of data first, and each row thereafter is the row north of the previous row. Most GIS applications write the array in a north-to-south fashion. The hydro model grid (Fullsom_hires.nc) is written north-to-south in accordance with GIS applications.

I looked over your python script. Although I did not test your code, but it does look like you are correctly reading, flipping, and applying the CRS to netCDF variables using xarray/rioxarray. 

There seems to have been a change at some point that breaks the ability of GIS applications such as QGIS to directly interpret the spatial metadata of the WRF-Hydro output files. I will have to test this a bit further because this is likely a small metadata change. If you wish to export a gridded WRF-Hydro variable to GeoTiff, you can use GDAL, which correctly interprets the projection and geotransform information. Here is an example of how you can export a GeoTiff on the command-line using GDAL:

gdal_translate NETCDF:202002010300.LDASOUT_DOMAIN1:ACCET ACCET.tif

 

I also looked at the frxst_pts in your Fulldom_hires.nc file and I see that the locations of the forecast points on the hydro grid are not representative of the points in the .csv file you provided. Are you certain that you are providing the correct CSV file to the script? I plotted the CSV points on this map (below) and I see where you are trying to locate these points. However, the placements on the grid seem to be wrong. Double check the CSV file to make sure you have the correct locations. Also, it looks like you the `-b False` is being interpreted incorrectly, so you can remove this from your command-line inputs and it should extend your flowline network downstream to your 4th forecast point.

Thanks,

Kevin

To view this discussion visit https://groups.google.com/a/ucar.edu/d/msgid/wrf-hydro_users/CA%2BSiuXKRpCcbqbdK_mMey6cH34DB%3DBEtiDt7LQbF8AZCkcxGSQ%40mail.gmail.com.

[Martin Nguyen]

Hi Kevin,

Thank you so much for your instructions. They are really helpful for me. I also checked again the forecast points and you were right. I changed them, checked again, and figured out that there are some other parameters in the hydro.namelist I think I should also changed accordingly: 

- DXRT = (dx or dy in namelist.wps) / (regridding factor in Build_Routing_Stack.py). In my case: DXRT = 200/4 = 50

- AGGFACTRT = regridding factor in Build_Routing_Stack.py. In my case: AGGFACTRT = 4

The frxst_pts in the Fulldom_hire.nc does not show the forecast points (like in the frxst_pts when I use  Examine_Outputs_of_GIS_Preprocessor.py. However, when I changed the above parameters and kept using that Fulldom_hires.nc, it still worked! Could you please help me check if my understanding here is correct? Also, for the FORCING data, should I also flip and apply crs (the same as the TOPOGRAPHY in the Fulldom_hires.nc) or I just need to keep it as it is?

I have updated the vers_003_streamflow.

Kind regards

Martin Nguyen

[Kevin Sampson]

Martin,

I am not sure I understand the issue with your forecast points. Did you zoom into the area containing the points? It can be hard to see one cell on a 50m grid when you are zoomed to the full domain. Your forecast points are clustered together near the basin outlet, so unless you apply a symbology with unique values, and zoom into the area, they can be hard to see.

With your forcing data, please keep it as-is. The model knows how to interpret the orientation of the gridded data, and has an understanding of the coordinate system from parameters stored in geogrid and wrfinput. So you do not need to make any changes to the forcing data for model execution, as long as the grids remain unchanged between your WRF simulation and your WRF-Hydro simulation. Adding the CRS and changing the orientation of the data is only needed for your visualization purposes. 

Thanks,

Kevin

[Martin Nguyen]

Hi Kevin,

Thank you so much for your information. I have tried some of more scenarios and they all work fine so far according to your instructions. About the issue with the forecast points, the frxst_pts in Fulldom_hires did not show the forecast point, but the Examine_Outputs_of_GIS_Preprocessor.py and the WRF-Hydro all show that the forecast points are included in the Fulldom_hires, so I guess that is fine.

I guess all my questions about the CHRTOUT_GRID1 of WRF-Hydro have already been solved. Just one more question about the WRF-Hydro results relating to flooding to do the analysis. Do you have any ideas about other results rather than streamflow that can help me show the flood information or to do the analysis? As for other hydrodynamic models like LISFLOOD-FP, they have the water surface elevation, and I tried to find it in the output lists here, but only the LAKE output has it.

Kind regards

Martin Nguyen

[Aubrey Dugger]

Hi Martin:

For your last question, the sfcheadsubrt field in RTOUT will have local inundation. WRF-Hydro does not currently do overbank river flooding - this needs to be done as a post-process (e.g., HAND methods).

Hope that helps.

Thanks!

Aubrey

[Martin Nguyen]

Hi Aubrey,

So sorry for the late reply and thank you so much for your last answer. I have tried the RTOUT and will have a look on HAND methods. I just have one more question. Is there any way that I can generate the CHANNELGRID only before doing Build_Routing_Stack.py? As if I can generate CHANNEL_GRID before, I can decide the forecast points easier to generate the forecast points csv file.

Kind regards

Martin Nguyen

[Kevin Sampson]

Martin,

Getting station points is an iterative process. If you keep your DEM and channel thresholds constant, then you can re-use the channelgrid from a previous run of Build_Routing_Stack.py to move your forecast point locations onto those grid cells. You can reference the GIS Pre-processing tools documentation, Section 9.1, for more information on this process.

Kevin

[Martin Nguyen]

Hi Kevin,

Thank you very much for your reply. I think that covers all the questions I have for now. If I have any further questions, I'll create a separate message.

Once, thank you to you, Aubrey, and Yongxin for all your help and patience. I really appreciate it.

Kind regards,

Martin Nguyen