# How to map emission inventory from lat&lon corrdinate to WRF model grid

Here, I have a regional emission inventory, which mainly contains the varibales of lat,lon,emission,and the emissionshape is (lat,lon).

For example :

1. ODIAC, the spatial resolution is 1*1 deg, and emission shape is (180,360)，units : gC/m2/d;

2. Hestia, the spatial resolution is 0.002*0.002 deg, and the emission shape is (521,560),units: Kg;

For now, I have to map the emission from the geographic coordinate (lat,lon) to WRF model grid, which may look like 99*99 with grid spacing of 3km,or 87*87 with grid spacing of 1km.

It seems that the process of mapping involves the calculation of area of old grid of emission inventory.

I am confused about how to do. Is there any efficient method to deal with it?

Thanks for all the help.

• The unit of the ODIAC emission inventory is given in gC/m2/d. You can interpolate this data set without considering the size of the grid cells. I do such remapping tasks with cdo. If you want to do the same for the Hestia data set with cdo, maybe the option setgridarea helps -- but I am not sure. But generally: yes, someone (the program or you) has to calculate the grid cell area. – daniel.heydebreck Mar 22 '18 at 14:28
• @daniel.neumann, thanks for your reply, I would try to learn the usage of cdo first. – Allen Zhang Mar 22 '18 at 15:21
• I might recommend creating a python script using the Basemap module (matplotlib.org/basemap) to convert both model and emissions geography into x and y coordinates, then interpolate the data from there. – BarocliniCplusplus Mar 22 '18 at 21:26
• @BarocliniCplusplus Is there any more detail example for it? – Allen Zhang Mar 23 '18 at 3:15
• Basemap is used to create maps. You can make a map identical to the WRF grid you use. Examples on how Basemap is commonly used can be found here: matplotlib.org/basemap/users/examples.html. I'd then recommend using the scipy interp module(docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html) to interpolate your ungridded data to the grid. – BarocliniCplusplus Mar 23 '18 at 14:43

I would suggest using cdo for your purpose. At least for variables, which values are independent of the grid cell size, one can use cdo.

# cdo

cdo (Climate Data Operators) is a command line program to process netCDF and GRIB files. It is developed an maintained by Uwe Schulzweida and colleagues at the Max-Planck Institute for Meteorology in Hamburg, Germany.

cdo provides different operators to process your data. You can chain several operators in one command. You always start with the call of cdo. Examples:

cdo OPERATOR1 -OPERATOR2 -OPERATOR3 inputfile1.nc inputfile2.nc outputfile.nc
cdo OPERATOR1 -OPERATOR2 inputfile1.nc -OPERATOR3 inputfile2.nc outputfile.nc


The first operator is just appended to cdo. All further appended operators need to have - as suffix.

## example

Assume, we have two files (file1.nc and file2.nc) with some ocean data (sst, salinity, u, v, ...). We want to calculate the difference between the sst of both files. This is done via

cdo diff -selname,sst file1.nc -selname,sst file2.nc outfile.nc


We could do this also in four steps:

cdo selname,sst file1.nc sst_of_file1.nc
cdo selname,sst file2.nc sst_of_file2.nc
cdo diff sst_of_file1.nc sst_of_file2.nc outfile.nc


# Interpolate data with cdo without considering grid cell size

cdo provides different operators for interpolation. REMAPBIL for bilinear interpolat, REMAPNN for nearest-neighbor interpolation, and so on. They are listed in the manual. I use remapbil here as an example.

cdo remapbil,NEW_GRID_DEF file_old_grid.nc file_new_grid.nc


The file_old_grid.nc is your input file and file_new_grid.nc is your output file. Here you find an example from my work at GitHub:

The grid definition on your input file file_old_grid.nc needs to be recognized by cdo. To test, whether it is recognized you can try

cdo griddes file_old_grid.nc


If it prints out the correct grid definition, then everything is fine. If not, you need to correct the definition of your lon-lat grid. Please have a look in the respective section below for details.

The NEW_GRID_DEF is a file, in which the target grid is described. If your target grid is a nicely defined grid, you can write this file by hand. E.g. it could look like (copied from the HTML manual of cdo):

gridtype = lonlat
xsize    =   60
ysize    =   30
xfirst   = -177
xinc     =    6
yfirst   =  -87
yinc     =    6


Please look for "CDO grid" and "grid description" in the cdo HTML Manual to fine more examples and explanation on it.

If you don't have a nicely defined grid, you can generate your NEW_GRID_DEF. This is, what I have often done, because our chemistry transport model setup needed a non-nicely defined grid.

cdo griddes file_with_wrf_grid_definition.nc > NEW_GRID_DEF


In file_with_wrf_grid_definition.nc, the model grid has to be defined and it needs to be recognized by cdo. Here you find an example from my work at GitHub:

## How to make cdo recognize

cdo assumes that your file conforms with the CF Conventions (Climate and Forecast Conventions). Please have a look in the respective chapter of the CF Conventions Document to get a general description.

A working grid definition could look as follows (as ncdump output):

dimensions:
time = UNLIMITED ; // (1 currently)
lay = 1 ;
lat = 78 ;
lon = 62 ;
variables:
float lat(lat, lon) ;
lat:long_name = "latitudes       " ;
lat:units = "degree north    " ;
lat:var_desc = "latitude (south negative)                                                       " ;
lat:cell_methods = "TSTEP, LAY: mean" ;
float lon(lat, lon) ;
lon:long_name = "longitudes      " ;
lon:units = "degree east     " ;
lon:var_desc = "longitude (west negative)                                                       " ;
lon:cell_methods = "TSTEP, LAY: mean" ;
float DATA_VAR(time, lay, lat, lon) ;
DATA_VAR:long_name = "some long var name" ;
DATA_VAR:units = "some unit" ;
DATA_VAR:var_desc = "some description" ;
DATA_VAR:coordinates = "lat lon" ;


The important aspect is that the data variable has an attribute coordinates, which value contains the two lon and lat dimensional variables (I attached a second example so that it becomes clear). Additionally, the standard_name and units attributes should be properly set at the dimensional variables.

dimensions:
time = UNLIMITED ; // (1 currently)
lay = 1 ;
lat = 78 ;
lon = 62 ;
variables:
float lat_1(lat, lon) ;
lat_1:long_name = "latitudes       " ;
lat_1:standard_name = "latitudes       " ;
lat_1:units = "degree north    " ;
lat_1:var_desc = "latitude (south negative)                                                       " ;
float lon_1(lat, lon) ;
lon_1:long_name = "longitudes      " ;
lon_1:standard_name = "longitudes      " ;
lon_1:units = "degree east     " ;
lon_1:var_desc = "longitude (west negative)                                                       " ;
float DATA_VAR(time, lay, lat, lon) ;
DATA_VAR:long_name = "some long var name" ;
DATA_VAR:standard_name = "the standard name" ;
DATA_VAR:units = "some unit" ;
DATA_VAR:var_desc = "some description" ;
DATA_VAR:coordinates = "lat_1 lon_1" ;


# Interpolate data with cdo with considering grid cell size

No idea how it works with cdo

• Hello,daniel. When thinking again this problem, since the unit of Hestia is Kg, I remapbil it to a 87*87 domain directly, it seems the result is a little strange. The maximum decreased from 169857 to 2771 Kg. Is it because of without considering the grid cell area or something else? However, when I first use gridboxsum,5,5 commond and then remapbil, the maximum of data just decrease to 124262 Kg. I am confused about that? – Allen Zhang Jun 11 '18 at 14:04
• On the other hand, It looks like we can use gridarea to calculate the cell_area, and then merge it to the file, then use exprwithnew_variable=original_variable/cell_area to get a new_variable with Kg/m2/h.But the shape of original_variable is (time,lat,lon),and the shape of cell_area is (lat,lon), do you have some experience about how to deal with it? – Allen Zhang Jun 11 '18 at 14:15
• it seems that use cdo div infile.nc -gridarea infile.nc outfile.nc can get the value of wanted new_variable, but I am not sure whether the whole process is right? since after remapbil, the maximum value would decrease from 4.46 to 0.073. – Allen Zhang Jun 12 '18 at 7:12

using R:

Alright,

As you have data with different resolution, I would suggest you to convert your data to raster and resample to match the grid cell of your wrf_inputs. Take care of mass conservation. Then convert your raster to a matrix, data.frame, or spatial feature of 'POLYGON' your emissions Then create wrf_chemi files

steps:

3. resample emissions
4. create data-frames of emissions
6. put your pollutants in wrf_chemis

Please, read the documentation of each R Package: raster, ratmos, sf, vein, eixport.

You can raster your data by:

library(raster)
devtools::install_github("ibarraespinosa/ratmos")
library(ratmos)
library(sf)
library(vein)
library(eixport)
r <- raster("emissions.tif") #or .nc, raster will use gdal to know
ri <- raster_wrf("wrfinput_d02")
spplot(ri, scales = list(draw = T)) # check plot
# Then it should be something like this
new_r <- crop(r, ri) # crop your data to the new raster
new_r <- resample(new_r, ri) #check options
df <- st_as_sf(rasterToPolygons(new_r))
# To create Array. Checkif you need to rotate your array
a <- GriddedEmissionsArray(df cols = ncol(ri), rows = nrow(ri), rotate = T)