爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 7453|回复: 7

[求助] python画图如何处理经纬度的错位(cmip6海温数据—nc文件)

[复制链接]

新浪微博达人勋

发表于 2022-5-4 11:06:06 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 wuliqin 于 2022-5-5 13:53 编辑

请教各位会使用python处理数据和画图的前辈、同伴们,当数据的结构是以下的类型,导致画出来的图经纬度错位了(经度错位90°,纬度从0开始),应该使用什么方法调整呢?我网上一直找不到解决方法,不知道有没有小伙伴遇到过跟我一样的问题然后有方法解决的,感谢指点!!

import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import xarray as xr
from matplotlib.colors import ListedColormap

file = 'F:/thetao_Omon_ACCESS-CM2_historical_r1i1p1f1_gn_201001-201412.nc'
a1 = xr.open_dataset(file)
print(a1)   # thetao  (time, lev, j, i) float32 ...

lat = (a1['j'][:])
lon = (a1['i'][:])
time = (a1['time'][:])
lev = (a1['lev'][:])
thetao = a1['thetao'][:][:, :, :, :]
print(thetao)
print(thetao.shape)   # (60, 50, 300, 360)

map = Basemap(lat_0=0, lon_0=0, llcrnrlat=-90., llcrnrlon=0., urcrnrlat=90., urcrnrlon=360., resolution='l')
lon, lat = np.meshgrid(lon, lat)    # np.meshgrid生成网格点坐标矩阵
x, y = map(lon, lat)                  # 将网格点坐标对应转换成地理经纬度

cmap = plt.cm.ocean_r        # 和以下这3行是截取想要的cmap
newcolors =cmap(np.linspace(0, 1, 256))
newcmap = ListedColormap(newcolors[10:190])

map.drawmeridians(np.arange(0, 360.001, 20), labels=[1, 0, 0, 1], linewidth=0)
map.drawparallels(np.arange(-90, 90.001, 10), labels=[1, 0, 0, 1], linewidth=0)

map.drawcoastlines(linewidth=0.2, color='black')           # 绘制海岸线
map.contourf(x, y, thetao[2, 5, :, :], levels=np.arange(0, 30, 2), cmap=newcmap, alpha=1.0, extend='both')    # 填色绘图
cbar = map.colorbar(location='bottom', pad="10%",)                   # 添加色标
cbar.set_label('℃')  # 添加色标标签

plt.show()
使用map.contourf(lon, lat, thetao[2, 5, :, :], levels=np.arange(-5, 35, 2), cmap=newcmap, extend='both', latlon=True)可以纠正经纬度的错位,但是出的图有问题

QQ图片20220504105341.png
QQ图片20220504105357.jpg
QQ图片20220505133720.png
QQ图片20220505133732.png
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2022-5-5 09:56:29 | 显示全部楼层
lat=a1.latitude
lon=a1.longitude
map.contourf(lon,lat,a1.sel(time=timestamp,level=level),latlon=True)
试试
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-5-5 13:48:21 | 显示全部楼层
气象小白01 发表于 2022-5-5 09:56
lat=a1.latitude
lon=a1.longitude
map.contourf(lon,lat,a1.sel(time=timestamp,level=level),latlon=Tr ...

太感谢您了,这样是可以的,只是第三句level出错了“keyerror”,我稍微改了一下就可以了。但是我不知道为什么,出的图在65°N以北特别奇怪,我用可视化软件panoply看了一下,数据没有问题,不知道哪一步出错了{:5_243:}
map.contourf(lon, lat, thetao[2, 5, :, :], levels=np.arange(-5, 35, 2), cmap=newcmap, extend='both', latlon=True)
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-5-6 09:11:16 | 显示全部楼层
不需要lon, lat = np.meshgrid(lon, lat) 和x, y = map(lon, lat) 。本身lat,lon就是二维网格了
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-5-6 09:14:03 | 显示全部楼层
用下面两句确定lat,lon
lat=a1.latitude
lon=a1.longitude
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-5-6 10:23:37 | 显示全部楼层
气象小白01 发表于 2022-5-6 09:14
用下面两句确定lat,lon
lat=a1.latitude
lon=a1.longitude

是的哦,是用的这个,代码如下,可是图在65°N以北还是有问题。

import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import xarray as xr
from matplotlib.colors import ListedColormap

file1 = 'F:/thetao_Omon_ACCESS-CM2_historical_r1i1p1f1_gn_201001-201412.nc'
a1 = xr.open_dataset(file1)
print(a1)   # thetao  (time, lev, j, i) float32 ...

lat = a1.latitude
lon = a1.longitude

time = (a1['time'][:])
lev = (a1['lev'][:])
thetao = a1['thetao'][:][:, :, :, :]
print(thetao)
print(thetao.shape)   # (60, 50, 300, 360)

map = Basemap(lat_0=0, lon_0=180, llcrnrlat=-90., llcrnrlon=0., urcrnrlat=90., urcrnrlon=360., resolution='l')

cmap = plt.cm.RdBu_r        # 和以下这3行是截取想要的cmap
newcolors =cmap(np.linspace(0, 1, 256))
newcmap = ListedColormap(newcolors[20:256])

map.drawmeridians(np.arange(-180, 180.001, 20), labels=[1, 0, 0, 1], linewidth=0)
map.drawparallels(np.arange(-90, 90.001, 10), labels=[1, 0, 0, 1], linewidth=0)
map.drawcoastlines(linewidth=0.2, color='black')           # 绘制海岸线
map.contourf(lon, lat, thetao[2, 5, :, :], levels=np.arange(-5, 35, 2), cmap=newcmap, extend='both', latlon=True)
cbar = map.colorbar(location='bottom', pad="10%",)                   # 添加色标
cbar.set_label('℃')  # 添加色标标签

plt.show()
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2022-5-6 20:31:11 | 显示全部楼层
thetao=a1.thetao
map.contourf(lon,lat,thetao.sel(time=?,lev=?))
试试这样?
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2022-5-7 08:47:07 | 显示全部楼层
气象小白01 发表于 2022-5-6 20:31
thetao=a1.thetao
map.contourf(lon,lat,thetao.sel(time=?,lev=?))
试试这样?

不行呢,有两个keyerror
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表