How does cross_section recognize projections? #1354
Comments
|
To answer your direct question:
As of MetPy 0.12, it does so by having a In your example, this was one of the issues you encountered, since MetPy won't recognize a Cartopy CRS instead of a CFProjection in the 'crs' coordinate. Another issue (as seen in your second screenshot) is that xarray 0.15.1 broke MetPy 0.12.0's accessor functionality. This was fixed in the recent MetPy 0.12.1 release, which I would recommend upgrading to. To resolve your issue now, I would recommend taking advantage of the convenience functionality in xarray and MetPy, and modifying your code to something like import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from datetime import datetime
from metpy.interpolate import cross_section
dtselect = datetime(2018, 7, 2, 18)
access_url = dtselect.strftime('https://psl.noaa.gov/thredds/dodsC/Datasets/NARR/pressure/air.%Y%m.nc')
ttds = xr.open_dataset(access_url).metpy.parse_cf()
ttx = ttds.sel(time=dtselect)
# Create a line going through NYC
start = (41.26, -73.95)
end = (40.04, -73.70)
cross = cross_section(ttx, start, end, steps=20).set_coords(('lon', 'lat'))
print(cross)(given that NARR is not a high-resolution dataset, it's probably a good idea to reduce the number of cross section steps from the default of 100 for the closely spaced start and end points.) |
jthielen
added
Area: Cross-sections
|
Thank you @jthielen for your response. I noticed that a new release of MetPy was released and it did finally fix the issues I was having with the MetPy example of cross_section.py script. However when I tried your fix for my specific dataset, I came across an issue at the I could be wrong, but this looks like I am using conda environments and I am attaching the environmental yml file I am using. Note that MetPy is being installed using pip since the latest version was not on conda-forge at the time I wrote it. For completeness here is a list of the versions of the packages being used for this script (I did Also thank you for your recommendation on the step size for cross_section. Once I get it working, I will definitely use less than 100 points, since it makes sense to do so for this resolution. |
|
I just tried running your recommended script on an updated xarray (0.15.1), and the same problem remains. Is there a way for me to use something like my original method, where I have an array of values, and then read only the metadata using |
|
Though, it still should be possible to work around it. I would think that, assuming in your original example, you can get import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from datetime import datetime
from siphon.catalog import TDSCatalog
from netCDF4 import Dataset
from metpy.interpolate import cross_section
base_url = 'https://psl.noaa.gov/thredds/catalog/Datasets/NARR/pressure/'
dtselect = datetime(2018, 7, 2, 18)
ttdscat = TDSCatalog(f'{base_url}/catalog.xml').datasets[f'air.{dtselect:%Y}{dtselect:%m}.nc'].access_urls['OPENDAP']
ttds = Dataset(ttdscat)
ttds = xr.open_dataset(xr.backends.NetCDF4DataStore(ttds)).metpy.parse_cf()
ttx = ttds.sel(time=dtselect)
# Create a line going through NYC
start = (41.26, -73.95)
end = (40.04, -73.70)
cross = cross_section(ttx, start, end, steps=20).set_coords(('lon', 'lat'))If that doesn't work, and it really is a memory issue so that you have to perform that index-based sub-setting with the NetCDF4 dataset...then it's a bit more involved workaround. It would then look something like this: import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from datetime import datetime
from siphon.catalog import TDSCatalog
from netCDF4 import Dataset, num2date, date2index, date2num
from metpy.interpolate import cross_section
from metpy.plots.mapping import CFProjection
base_url = 'https://psl.noaa.gov/thredds/catalog/Datasets/NARR/pressure/'
dtselect = datetime(2018, 7, 2, 18)
ttdscat = TDSCatalog(f'{base_url}/catalog.xml').datasets[f'air.{dtselect:%Y}{dtselect:%m}.nc'].access_urls['OPENDAP']
ttds = Dataset(ttdscat)
# Now we select our time of interest from the data for the temeprature variable.
dtunits = ttds.variables['time'].units # time epoch
t0time = date2num(dates = dtselect, units = dtunits)
t0index = int(np.where(ttds.variables['time'][:] == t0time)[0]) # start time index
times = num2date(ttds.variables['time'][t0index], units = dtunits) # Datetimes, from file
tt = ttds.variables['air'][t0index, :, :, :]
# Extract the level of interest and the coordinates
lat = ttds.variables['lat'][:]
lon = ttds.variables['lon'][:]
level = ttds.variables['level'][:]
dataproj = CFProjection({
key: ttds.variables['Lambert_Conformal'].getncattr(key) for key in ttds.variables['Lambert_Conformal'].ncattrs()
})
ttx = xr.Dataset({'temperature': (['level', 'y', 'x'], tt),
'lat': (['y', 'x'], lat),
'lon': (['y', 'x'], lon)},
coords = {'x': ttds.variables['x'][:],
'y': ttds.variables['y'][:],
'time': dtselect,
'level': level,
'crs': dataproj})
ttx = ttx.metpy.parse_cf()
# Create a line going through NYC
start = (41.26, -73.95)
end = (40.04, -73.70)
cross = cross_section(ttx, start, end).set_coords(('lon', 'lat')) |
jthielen
added
Status: Upstream
|
So the problem is that when we do Luckily, we can work around this by changing the order of some operations: import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from datetime import datetime
from metpy.interpolate import cross_section
dtselect = datetime(2018, 7, 2, 18)
access_url = dtselect.strftime('https://psl.noaa.gov/thredds/dodsC/Datasets/NARR/pressure/air.%Y%m.nc')
ttds = xr.open_dataset(access_url)
subset = ttds.sel(time=dtselect)
ttx = subset.metpy.parse_cf()
# Create a line going through NYC
start = (41.26, -73.95)
end = (40.04, -73.70)
cross = cross_section(ttx, start, end, steps=20).set_coords(('lon', 'lat'))
print(cross)That works for me (and works quite quickly!). @jthielen any idea if there's anything we can do to make the internal call to |
|
@dopplershift Good catch, and thanks for the correction! I tripped myself up with the different versions at play. I realized now that my initial and second suggestion only worked on the MetPy 0.12.x without But, this subtle point raises a big issue. Before the problematic call to either ttds['air'].variable._datawhich is definitely not something that Pint can handle wrapping at the moment, and even if it was, xarray wouldn't be able to handle it being inside an external wrapper right now either. I'm now at a bit of a loss though about the best way forward on this would be, however, since this is now a somewhat major downside to forcing wrapping in Quantity that I neglected to consider in #1304 / #1325. To at least put some ideas out there, should we:
|
|
Also, side note, @dopplershift I realize I've been a bit too overeager with the labels, leading to some inaccuracies that have needed to be corrected later. I'll dial it back a bit. |
|
After thinking about this a bit more, I think the best way forward is to:
Does this seem reasonable? |
|
That seems reasonable to me. The only question I'd double check with is: does any of this interfere with a dtype units (NEP-41) world? I know were low on details as of yet, but I just want to make sure we're not going down any path that obviously interferes. |
|
For documenting purposes, here is the the play-by-play on how each code performed and my takeaways at the end. Thank you @jthielen and @dopplershift for your feedback!! @jthielen's SuggestionsSuggestion 1: "Simple"
Result: Same error occurred. Suggestion 2: More Involved
Result: This code worked!! I added a line to print @dopplershift's SuggestionSuggestion 3: Change Order of OperationsThat clarification on
Result: This also worked for me. The output: My TakeawaysThank you @jthielen and @dopplershift! I have a way forward now, and look forward to producing the corresponding plots for these cross-sections (similar to the one your MetPy cross_section.py example). My first takeaway is that I can be careful about the order of my operations so that My other takeaway is that projection information can be supplied in an 'ad hoc' manner using the If my WRF is running a "Lambert Conformal" projection, then I would just use the identical keys here and swap the values. But if I was using a different projection I would probably have to find the correct parameters for that projection. After some digging, I found that those keys are found, verbatim in CF Metdata in Table 4, under Lambert Conformal (LCC-1SP). The code snippets y'all provided, with this link for CF Metadata, seems to be the complete answer to my question. My apologies if this link (or very similar material) was already referenced somewhere in the code, and I am just ignorant about it. |
I would guess that it does not interfere implementation-wise (or at least change anything with respect to that). I would think altering the xarray accessor (what this comes down to) to operate on a dtype-backed-unit array instead of a Pint Quantity as its preferred unit array type would be almost trivial compared to the modifications needed throughout the rest of the library. Where it would be good to keep this in mind, however, is in writing documentation, making sure to not give impressions that this way of working with units and the workarounds still needed are how it will be from now on. |
|
Also, thank you @thedmv for the great synopsis! I would have like to have said in regards to your WRF takeaway that there are helpers in place for that already (that was my plan earlier this year: https://ams.confex.com/ams/2020Annual/meetingapp.cgi/Paper/369355), but they aren't quite there yet. For now, I'd just want to point out https://gist.github.com/jthielen/8881d32d08d625c75c906a7e8ad7583f and #1089, which may give you a head start. |
Problem Description
I am attempting to use the cross_section function for NARR data that I have downloaded from PSL @ NOAA. I keep getting errors around how handle coordinates properly.
My Code
My error
Looking at cross_section.py for help
I was hoping to learn from the example in order to fix my problems, but I can't get past this error from the example script cross_section.py.
I'm not sure where to start to begin to address my problem with defining the CRS for my data, and I was hoping to seek guidance from your cross_section.py code, but now I'm stuck with nothing working for me.
-David
Python version: Python 3.7.6
MetPy version: 0.12.0
I am using Python on Windows using Miniconda 3.
The text was updated successfully, but these errors were encountered: