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? Specifically http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Basics/GEOGRID/geogrid_namelist.htm and http://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 to i_parent_start and j_parent_start? How is this related to dx and dy and ref_lat and ref_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) Run geogrid.exe to generate the parent grid. Since you don't know the exact location of the nest yet, set max_dom = 1 in namelist.wps. This will generate a file called geo_em.d01.nc.

2) Look at the geo_em.d01.nc file to find the right indices for your child domain. i_parent_start and j_parent_start refer to the x and y indices on the parent grid at which the southwest corner of the child grid will be positioned. XLONG_M and XLAT_M are 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)

where xlong_m and xlat_m are 2-dimensional arrays that you read from the grid, and lat0 and lon0 are 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 edit namelist.wps again, set the i_parent_start and j_parent_start to the values that you calculated in step 3, set max_dom = 2, and re-run geogrid.exe. The child domain file geo_em.d02.nc should be created.

4) Look at the geo_em.d02.nc file. Repeat the procedure until happy with the domain location.

About the parent_grid_ratio parameter. 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 for parent_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 = 3 is 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 calculates i_parent_start and j_parent_start. It is in Java and as people are aware Java does not have functions such as minloc or Python's unravel_index. So I divided my program into three steps.

1. Read in the xlat_m and xlong_m arrays using java netcdf.

2. Subtracting desired lat and desired lon from each value of the 2-dimensional xlat_m and xlong_m arrays and squaring the difference.

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

There is a tool that comes with WRF called the WRF Domain 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 (either geo_em or wrfout NetCDF 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)