Prepare the emission inventory file for MPAS-CMAQ

MPAS-CMAQ
,
1

Hello!

I am trying to interpolate the EDGAR v8.1 emission inventory onto a 120 km uniform grid of MPAS. I used the following code to extract the MPAS grid information:

def create_mpas_grid(mpas_file,proj,filename):
    ds = xr.open_dataset(mpas_file)

    indexToCellID = ds['indexToCellID'].values

    verticesOnCell = ds['verticesOnCell'].values
    nEdgesOnCell = ds['nEdgesOnCell'].values

    lonVertex = np.degrees(ds['lonVertex'].values)
    latVertex = np.degrees(ds['latVertex'].values)

    polygons = []
    for i in range(len(indexToCellID)):
        vertex_indices = verticesOnCell[i, :nEdgesOnCell[i]] - 1

        poly_lons = lonVertex[vertex_indices]
        poly_lats = latVertex[vertex_indices]
        polygons.append(Polygon(zip(poly_lons, poly_lats)))

    mpas_grid = gpd.GeoDataFrame({'cellID': indexToCellID}, geometry=polygons)
    mpas_grid.crs = proj
    if os.path.exists(filename+".shp"): os.remove(filename+".shp")
    if os.path.exists(filename+".shx"): os.remove(filename+".shx")
    if os.path.exists(filename+".dbf"): os.remove(filename+".dbf")
    if os.path.exists(filename+".prj"): os.remove(filename+".prj")
    if os.path.exists(filename+".qpj"): os.remove(filename+".qpj")
    if os.path.exists(filename+".cpg"): os.remove(filename+".cpg")
    mpas_grid.to_file(filename+".shp")

    return mpas_grid

# mpas grid
if os.path.exists('mpasgrid.shp'):
    grid2 = gpd.read_file('mpasgrid.shp')
else:
    mpas_file = 'x1.40962.grid.nc'
    proj = 'epsg:4326'
    grid2 = create_mpas_grid(mpas_file,proj,'mpasgrid')

Then, I checked the coverage of the MPAS grid generated from the extracted information, and found that it is not fully closed globally.

print(grid2.total_bounds)

图片

I would like to know if this is an issue with the MPAS grid itself or if my method of extracting the grid information is too rough? If it is the latter, could you teach me a better way to do it?
I look forward to your reply! Thank you!

Sincerely,
Pan.

2

The vertices of the bounding box are likely correct given the extents of the 120 mesh. While a polygon overlap weighting method works for regridding, there are readily available tools for interpolation from a lat-lon grid to MPAS mesh that may be more efficient for your needs. For example, the pyremap tool (pyremap — pyremap main documentation) can generate a regridding matrix between MPAS, projected, and lat-lon grids.

2 Likes
3

Hello @james.beidler

Your answer was very helpful to me! I will give it a try.But I would like to know, since emissions are different from variables like temperature and wind speed that can be interpolated using nearest neighbors and are strongly correlated with grid area, can this tool easily perform interpolation for such variables?

Wish you all the best,
Pan.

4

Yes, pyremap can interface with ESMF for re-gridding. The ESMF interpolation methods are described here and include conservative methods: 2 Command Line Tools

5

Hello,

Thank you for your explanation. It was very helpful to me!

Wish you all the best,
Pan.