Search code examples
pythonpython-xarraycartopymetpy

Question about the coordinates of metpy cross section


I followed the code to draw the cross section figure on this site. (https://unidata.github.io/MetPy/latest/examples/cross_section.html#sphx-glr-examples-cross-section-py)

In the example, running the following code will produce this result(cross). cross = cross_section(data, start, end).set_coords(('lat', 'lon')) enter image description here

However, unlike the example, the x,y coordinates of the result value after running the cross section code have changed. And the longitude and latitude are fixed. enter image description here I didn't understand why this result came out like this.

my code

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import numpy as np
import os
import xarray as xr

import metpy.calc as mpcalc
from metpy.cbook import get_test_data
from metpy.interpolate import cross_section

os.chdir ("D:/PPL/CrossSection/20190404/")
ncfile = r'./2019040400pres.nc'
NC = xr.open_dataset(ncfile)

pres_list = ['1000','975','950','925','900','875','850','800','750','700','650','600','550','500','450','400','350','300','250','200','150','100','70','50'];
p_list = np.array(pres_list, dtype = np.float64)

NC = NC.metpy.parse_cf().squeeze()

NC = NC.assign_coords(isobaric = p_list).set_coords('isobaric')

My code had to be analyzed based on the nc file divided by the isobaric layer, so I combined the values divided for each layer. like this,

tmp1 = NC['TMP_1000mb']
tmp2 = NC['TMP_975mb']
...
tmp23 = NC['TMP_70mb']
tmp24 = NC['TMP_50mb']

TMP = xr.concat([tmp1,tmp2],'isobaric')
TMP = xr.concat([TMP,tmp3],'isobaric')
...
TMP = xr.concat([TMP,tmp23],'isobaric')
TMP = xr.concat([TMP,tmp24],'isobaric')

TMP['isobaric']=p_list

NC['Temperature'] = TMP

So, i print NC. The result is

<xarray.Dataset>
Dimensions:            (isobaric: 24, x: 602, y: 781)
Coordinates:
    * y                  (y) float64 0.0 1.5e+03 3e+03 ... 1.168e+06 1.17e+06
    * x                  (x) float64 0.0 1.5e+03 3e+03 ... 9e+05 9.015e+05
    latitude           (y, x) float64 32.26 32.26 32.26 ... 42.94 42.94 42.93
    longitude          (y, x) float64 121.8 121.9 121.9 ... 132.5 132.5 132.5
    time               datetime64[ns] 2019-04-04
    metpy_crs          object Projection: latitude_longitude
    * isobaric           (isobaric) float64 1e+03 975.0 950.0 ... 100.0 70.0 50.0
Data variables:
    ...
    ...
    Temperature        (isobaric, y, x) float32 282.87515 282.87515 ... nan nan

Attributes:
    Conventions:          CF-1.0
    History:              created by wgrib2
    GRIB2_grid_template:  30

However, when I run cross section, x,y changes to start and end, and the original longitude and latitude are fixed.

start = (37.5, 126.63)
end = (37.8 ,128.86)
cross = cross_section(NC, start, end)
#I didn't write .set_codes here because unlike the example, coordinates are given 
#latitude and longitude.

print(cross)
<xarray.Dataset>
Dimensions:            (index: 100, isobaric: 24)
Coordinates:
    latitude           (index) float64 32.26 32.26 32.26 ... 32.26 32.26 32.26
    longitude          (index) float64 121.8 121.8 121.8 ... 121.8 121.8 121.8
    time               datetime64[ns] 2019-04-04
    metpy_crs          object Projection: latitude_longitude
    x                  (index) float64 126.6 126.7 126.7 ... 128.8 128.8 128.9
    y                  (index) float64 37.5 37.5 37.51 37.51 ... 37.79 37.8 37.8
  * index              (index) int32 0 1 2 3 4 5 6 7 ... 92 93 94 95 96 97 98 99
  * isobaric           (isobaric) float64 1e+03 975.0 950.0 ... 100.0 70.0 50.0
Data variables:
    ...
    ...
    Temperature        (isobaric, index) float64 282.9 282.9 ... 209.6 209.6
Attributes:
    Conventions:          CF-1.0
    History:              created by wgrib2
    GRIB2_grid_template:  30

I want x,y to be maintained and latitude and longitude to be in the start and end ranges as in the example. I don’t really understand why this is happening.

Thank you for your attention.

For reference, I attach the data source and initial data. Data source : LDAPS Data from KMA

Full NC data(initial data)

<xarray.Dataset>
Dimensions:            (time: 1, x: 602, y: 781)
Coordinates:
* y                  (y) float64 0.0 1.5e+03 3e+03 ... 1.168e+06 1.17e+06
* x                  (x) float64 0.0 1.5e+03 3e+03 ... 9e+05 9.015e+05
latitude           (y, x) float64 ...
longitude          (y, x) float64 ...
* time               (time) datetime64[ns] 2019-04-04
Data variables:
DZDT_1000mb        (time, y, x) float32 ...
DZDT_975mb         (time, y, x) float32 ...
...
DZDT_70mb          (time, y, x) float32 ...
DZDT_50mb          (time, y, x) float32 ...
UGRD_1000mb        (time, y, x) float32 ...
VGRD_1000mb        (time, y, x) float32 ...
UGRD_975mb         (time, y, x) float32 ...
VGRD_975mb         (time, y, x) float32 ...
...
UGRD_70mb          (time, y, x) float32 ...
VGRD_70mb          (time, y, x) float32 ...
UGRD_50mb          (time, y, x) float32 ...
VGRD_50mb          (time, y, x) float32 ...
HGT_1000mb         (time, y, x) float32 ...
HGT_975mb          (time, y, x) float32 ...
...
HGT_70mb           (time, y, x) float32 ...
HGT_50mb           (time, y, x) float32 ...
TMP_1000mb         (time, y, x) float32 ...
TMP_975mb          (time, y, x) float32 ...
...
TMP_70mb           (time, y, x) float32 ...
TMP_50mb           (time, y, x) float32 ...
var0_1_194_1000mb  (time, y, x) float32 ...
var0_1_194_975mb   (time, y, x) float32 ...
...
var0_1_194_70mb    (time, y, x) float32 ...
var0_1_194_50mb    (time, y, x) float32 ...
RH_1000mb          (time, y, x) float32 ...
RH_975mb           (time, y, x) float32 ...
...
RH_70mb            (time, y, x) float32 ...
RH_50mb            (time, y, x) float32 ...
Attributes:
Conventions:          CF-1.0
History:              created by wgrib2
GRIB2_grid_template:  30

Solution

  • Notice that the metpy_crs identified in your dataset from .metpy.parse_cf() is given as Projection: latitude_longitude, which is not correct, as your dataset has 2D latitude/longitude coordinates with 1D y/x coordinates in a projected grid space. This is likely triggered by your dataset's metadata not being CF-compliant. At it is now, because of this, MetPy is treating your x coordinate as if it were longitude and y as if it were latitude (notice the ranges in those coordinates), even though they are not. The outcome here is definitely far from ideal (no error message or anything), so I'd recommend submitting an issue on MetPy's issue tracker.

    To directly fix your issue, you will need to manually specify the CRS/projection of your data. To do so (in MetPy v1.0 and later), use the assign_crs method on MetPy's Dataset accessor instead of using parse_cf (which only works for CF-compliant datasets).