parse_cf does not work when multiple grid mappings are present

[djhoese]

This bug was mentioned on gitter when talking to @dopplershift. I am beginning more in-depth designing of the geoxarray library which will likely borrow a lot from metpy's work with CF Dataset/DataArray objects. One concern I had was handling NetCDF files with multiple grid_mapping variables. It turns out this is not possible when parse_cf is used without arguments in one of these cases. See also geoxarray/geoxarray#2

Here is an example:

import xarray as xr
import numpy as np
import metpy

x = xr.DataArray(np.arange(5), attrs={'units': 'm'}, dims=('x',))
y = xr.DataArray(np.arange(5), attrs={'units': 'm'}, dims=('y',))
x2 = xr.DataArray(np.arange(10), attrs={'units': 'm'}, dims=('x2',))
y2 = xr.DataArray(np.arange(10), attrs={'units': 'm'}, dims=('y2',))                                                                                                                     

ds = xr.Dataset({
    'ds1': xr.DataArray(np.zeros((5, 5)), dims=('y', 'x'), coords={'y': y, 'x': x}, attrs={'grid_mapping': 'gm1'}),
    'ds2': xr.DataArray(np.zeros((5, 5)), dims=('y', 'x'), coords={'y': y, 'x': x}, attrs={'grid_mapping': 'gm1'}), 
    'ds3': xr.DataArray(np.zeros((10, 10)), dims=('y2', 'x2'), coords={'y2': y2, 'x2': x2}, attrs={'grid_mapping': 'gm2'}),
    'gm1': xr.DataArray([], attrs={'grid_mapping_name': 'lambert_conformal_conic'}), 
    'gm2': xr.DataArray([], attrs={'grid_mapping_name': 'transverse_mercator'})})

new_ds = ds.metpy.parse_cf()

This produces:

~/anaconda/envs/polar2grid_py36/lib/python3.6/site-packages/metpy/xarray.py in parse_cf(self, varname, coordinates)
    156         # If no varname is given, parse the entire dataset
    157         if varname is None:
--> 158             return self._dataset.apply(lambda da: self.parse_cf(da.name,
    159                                                                 coordinates=coordinates))
    160 
...
~/anaconda/envs/polar2grid_py36/lib/python3.6/site-packages/xarray/core/merge.py in unique_variable(name, variables, compat)
     81                                  'objects to be combined:\n'
     82                                  'first value: %r\nsecond value: %r'
---> 83                                  % (name, out, var))
     84             if combine_method:
     85                 # TODO: add preservation of attrs into fillna

MergeError: conflicting values for variable 'crs' on objects to be combined:
first value: <xarray.Variable ()>
array(<metpy.plots.mapping.CFProjection object at 0x1a27d70160>, dtype=object)
second value: <xarray.Variable ()>
array(<metpy.plots.mapping.CFProjection object at 0x1a27cefba8>, dtype=object)

This is happening because each DataArray inside the Dataset is being parsed separately and each one gets a crs coordinate added. When xarray gets to the point that it merges the DataArray objects together in a Dataset it can't handle merging the crs coordinates with different values.

[jthielen]

I'd take the git blame for this. MetPy's parse_cf() started off only parsing out a single variable, so there was no need to worry about conflicts in the crs coordinate. The naive dataset parsing in #871 was added as a convenience, but this issue demonstrates that more careful thought is needed for fully parsing datasets.

What would be the intended behavior here? Have a unique crs1, crs2, ... for each different grid_mapping in the dataset with some internal logic (other than by name like it is currently) for the helpers (like da.metpy.cartopy_crs) to access the proper CFProjection? Or some other approach?