How to get cell indices from latitude and longitude in WRF model grids?

Question

I am using Weather Research and Forecasting model WRF and the WRF Preprocessing System WPS and could not get a satisfying answer from wrfhelp. Can anyone guide me in transforming latitude and longitude to grid points especially for nested grids? Specificallyhttp://www2.mmm.ucar.edu/wrf/OnLineTutorial/Basics/GEOGRID/geogrid_namelist.htmandhttp://www2.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm#_Description_of_the_1

How does one transform latitude and longitude toi_parent_startandj_parent_start? How is this related todxanddyandref_latandref_lon? A simple example and an explanation would be very helpful.

Answer 1

I've been a WRF user for almost 5 years now, and contributed code to a recent public release. I am not aware that WPS (WRF Preprocessing System) has such a tool that takes in the grid and point coordinates and returns the appropriate index. However, it is very straightforward to do so yourself. Some suggest using an external library, I think that may be an overkill for such a simple task. Here is what you need to do:

1) Rungeogrid.exeto generate the parent grid. Since you don't know the exact location of the nest yet, setmax_dom = 1innamelist.wps. This will generate a file calledgeo_em.d01.nc.

2) Look at thegeo_em.d01.ncfile to find the right indices for your child domain.i_parent_startandj_parent_startrefer to the x and y indices on the parent grid at which the southwest corner of the child grid will be positioned.XLONG_MandXLAT_Mare the longitude and latitude grids of the mass (pressure) points. Using a programming language of choice, find the grid cell that is closest to your desired location for the child nest corner. This is typically done by looking for the minimum value of distance between desired location and all the points on the grid. For example, in Fortran, you can do something like:

integer :: i_parent_start,j_parent_start integer,dimension(2) :: coords coords = minloc((lon0-xlong_m)**2+(lat0-xlat_m)**2) i_parent_start = coords(1) j_parent_start = coords(2)

wherexlong_mandxlat_mare 2-dimensional arrays that you read from the grid, andlat0andlon0are the desired coordinates of the child nest southwest corner. Similarly, if you use Python, you could do:

import numpy as np j_parent_start,i_parent_start = np.unravel_index(\ np.argmin((lon0-xlong_m)**2+(lat0-xlat_m)**2),xlon_m.shape) # Add one because WRF indices start from 1 i_parent_start += 1 j_parent_start += 1

3) Now editnamelist.wpsagain, set thei_parent_startandj_parent_startto the values that you calculated in step 3, setmax_dom = 2, and re-rungeogrid.exe. The child domain filegeo_em.d02.ncshould be created.

4) Look at thegeo_em.d02.ncfile. Repeat the procedure until happy with the domain location.

About theparent_grid_ratioparameter. This is an integer factor of child grid refinement relative to the parent grid. For example, if set to 3, and parent grid resolution is 12 km, the child grid resolution will be 4 km. Odd values forparent_grid_ratio(3, 5, etc.) are recommended because for even values, interpolation errors arise due to the nature of Arakawa C-grid staggering.parent_grid_ratio = 3is the most commonly used value, and recommended by myself.

Answer 2

This question might be better suited for GIS stack exchange... but I'll take a shot:

dx and dy are the distance per grid cell in the projected grid units, which is part of your grid definition. ref_lat and ref_lon are part of your projection definition. i_parent_start and j_parent_start are the horizontal grid cell numbers of the parent grid, where your nested grid sits. I believe WRF grids are oriented so that i,j of 0,0 are the bottom left of your grid (e.g. a typical x-y cartesian coordinate system) and the first pixel position is defined to be at the 1,1 coordinate.

In order to create a shapefile of the grid projection, I would personally use the Spatial Allocator Tool from UNC (https://www.cmascenter.org/sa-tools/). I believe the same tool will map points to your grid as well. However, if you want to map things to check the results, you could use your projection definition to build a custom transformation in a GIS program, project your lat/longs, and then map everything together.

Answer 3

I am posting my code that calculatesi_parent_startandj_parent_start. It is in Java and as people are aware Java does not have functions such asminlocor Python'sunravel_index. So I divided my program into three steps.

1. Read in thexlat_m and xlong_marrays using java netcdf.

2. Subtracting desired lat and desired lon from each value of the 2-dimensionalxlat_m and xlong_marrays and squaring the difference.

3. Then finding the indicesiandjof the smallest element of the 2-dimensional array (not as a flattened 1-d array).

有一个工具,伴随着WRF WRFDomain Wizard that creates the parent and child grids but the issue with this is that you can't see the topography and land use that is so essential for high resolution simulations. In order to view the topography I have another tool that converts the NetCDF file into a KML file and project that on Google Earth. The enclosed image of the child grid is also presented. As explained above for coarse grids this is not so much an issue as it is for nested grids where in the size of the grid itself may not be more than a few hundred meters.

double desiredLatitude = 27; double desiredLongitude = 77; double[][] xlatMArray = getXLatMArray(dataFile); double[][] xlongMArray = getXLongMArray(dataFile); int a = xlatMArray.length; double[][] array = new double[a][a]; for (int i =0 ; i < a;i++) { for (int j =0;j < a ;j++) { double latitudeDiff = desiredLatitude - xlatMArray[i][j]; double longitudeDiff = desiredLongitude - xlongMArray[i][j]; array[i][j]= Math.pow(latitudeDiff,2) + Math.pow(longitudeDiff,2); } } double min = array[0][0]; double max = array[0][0]; int[] minIndex = {0,0}; int[] maxIndex = {0,0}; for (int i =0; i < array.length ; i++) { for (int j =0 ; j < array.length ; j++) { if (array[i][j] < min) { min = array[i][j]; minIndex[0] = i; minIndex[1] = j; } if (array[i][j] > max) { max = array[i][j]; maxIndex[0] = i; maxIndex[1] = j; } } }

Answer 4

Try NCL. It has a package for WRF-related functions, including the bidirectional grid-coordinates conversion routine.

To convert coordinate pairs to grid indices, follow the documentation here:wrf_user_ll_to_ij. A backward (i.e. gird to coordinate) routine can be found here:wrf_user_ij_to_ll. They are quite handy.

For your scenario, you may load your parent domain files (eithergeo_emorwrfoutNetCDF files), and then specify the child domain starting coordinate, you will get the results.

Here are some examples from the official documentation:

Conversion from coordinate to grid indices:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" a = addfile("wrfout_d01_2000-01-24_12:00:00.nc","r") loc = wrf_user_ij_to_ll(a, 10., 30., True) print("lon location is: " + loc(0)) print("lat location is: " + loc(1))

Conversion from grid indices to coordinate:

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl" a = addfile("wrfout_d01_2000-01-24_12:00:00.nc","r") res = True res@returnInt = False ; return real values loc = wrf_user_ll_to_ij(a, -100.0, 40.0, res) print("X/Y location is: " + loc)