nc文件数据处理后如何提取经纬度坐标和某段时间进行作图

lynpatty

之前一直都在做数据处理，没有过相关可视化的经验，代码敲下来以后想要把经纬度一起输出进行画图可是想来想去都不对，想请大家看看代码。如何把对应MD的格点输出，画出41年里某一个月的图呢？
import os
import time
from netCDF4 import Dataset as nc
import numpy as np
import pandas as pd
from tqdm import tqdm
import matplotlib.pyplot as plt###引入库包
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker


files = os.listdir('/downloads/Tmax/CPC_tmax') # read data

# initialization
tmax = np.zeros((41, 400, 40, 40)).astype('float64')

# read data
index = 0
year_marks = [] # 用来记录一年有多少天（主要是有闰年）
for f in tqdm(files):
        if '.nc' not in f:
                continue
        temp = nc('/downloads/Tmax/CPC_tmax/{}'.format(f))
        temp_tmax = temp.variables['tmax'][:,210:250, 240:280]
        year_marks.append(temp_tmax.shape[0])
        tmax[index,15:15+temp_tmax.shape[0],:,:] = temp_tmax.copy()
        index = index + 1

for i in range(40):
        ylen = year_marks
        tmax[i, ylen+15:400,:,:] = tmax[i+1, 15:400-ylen,:,:] # 将后年的开头补到今年的结尾
        tmax[i+1, 0:15,:,:] = tmax[i, ylen:ylen+15,:,:] #将今年的结尾补到后年的开头
tmax[tmax<-0.1] = np.nan

# 计算热浪事件
# 先计算每个格点的90分位线
print("90 percentile calculation starting... ")
threshold90 = np.zeros((400,40,40))
for i in tqdm(range(40)):
        for j in range(40):
                for k in range(15, 15+366):
                        data = tmax[:,k-15:k+16,i,j].copy()
                        data = data[~np.isnan(data)] # 去掉缺省值
                        if data.shape[0] >=50:
                                threshold90[k,i,j] = np.percentile(data, 90, axis=0)
indexes = [i for i in range(40)]

# 再对每个格点判断有没有超过90分位线
tmax_flags = np.zeros((41,400,40,40))
for i in tqdm(range(40)):
        for j in range(40):
                for y in range(41):
                        for k in range(15, 15+year_marks[y]):
                                if threshold90[k,i,j] != 0:
                                        if tmax[y,k,i,j] > threshold90[k,i,j]:
                                                tmax_flags[y, k,i,j] = 1

for y in range(41): # 因为之后要滑动窗口找热浪事件，所以这里也要处理头尾补全
        if y != 0:
                tmax_flags[y,14,:,:] = tmax_flags[y-1,15+year_marks[y]-1,:,:] # 给1-41年补头
        if y != 40:
                tmax_flags[y,year_marks[y],:,:] = tmax_flags[y+1,15,:,:] # 给0-40年补尾

# 计算 threshold25 和 threshold75
print("Calculating threshold25 & threshold75...")
# 1*41
threshold25 = np.zeros((40,40))
threshold75 = np.zeros((40,40))
for i in range(40):
        for j in range(40):
                threshold_temp = []
                for y in range(41):
                        threshold_temp.append(np.nanmax(tmax[y,15:15+year_marks[y],i,j]))
                threshold_temp = np.array(threshold_temp)
                threshold25[i,j] = np.percentile(threshold_temp,25, axis=0)
                threshold75[i,j] = np.percentile(threshold_temp,75, axis=0)
print("Done.")

print("Calculating Heatwaves...")'
timeseries_dict = {}
for y in range(41):
        num_event = 0
        md_values = []
        output = open(/downloads/Tmax/step5/{}_md.txt'.format(1979+y), 'w')
        for i in range(40):
                for j in range(40):
                        zeros = []
                        diff = (threshold75[i,j] - threshold25[i,j])
                        for k in range(year_marks[y]):
                                if tmax_flags[y,k+15,i,j] == 0:
                                        zeros.append(k+15)
                        for l in range(len(zeros)):
                                if l == 0:
                                        continue
                                else:
                                        interval = zeros[l]-zeros[l-1]
                                        if interval > 3:
                                                num_event += 1
                                                md_value = 0
                                                for kk in range(zeros[l-1]+1, zeros[l]):
                                                        if tmax[y,kk,i,j] > threshold25[i,j]:
                                                                print("hello?")
                                                                md_value += (tmax[y,kk,i,j] - threshold25[i,j]) / diff
                                                md_values.append(md_value)
        for i in md_values:
                output.write("{}\n".format(i))
        output.close()
        timeseries_dict[1979+y] = num_event
print("Done.")
timeseries_pd = pd.DataFrame.from_dict(timeseries_dict, orient='index')
timeseries_pd.to_csv("step5.csv")


proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(4, 4), dpi=550)  # 创建画布
extent=[120,140,-35,-15]
ax = fig.subplots(1, 1, subplot_kw={'projection': proj})  # 创建子图
ax.set_extent(extent)
ax.add_feature(cfeature.LAND)####添加陆地######
ax.add_feature(cfeature.COASTLINE,lw=0.3)#####添加海岸线#########
ax.add_feature(cfeature.RIVERS,lw=0.25)#####添加河流######
ax.add_feature(cfeature.LAKES)######添加湖泊#####
ax.add_feature(cfeature.BORDERS, linestyle='-',lw=0.25)###
ax.add_feature(cfeature.OCEAN)######添加海洋########
ax.add_feature(cfeature.OCEAN.with_scale('10m'))
ax.add_feature(cfeature.LAND.with_scale('10m'))
ax.add_feature(cfeature.RIVERS.with_scale('10m'),lw=0.6)
ax.add_feature(cfeature.LAKES.with_scale('10m'))
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle='-',lw=0.6)
ax.add_feature(cfeature.COASTLINE.with_scale('10m'),lw=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.2, color='k', alpha=0.5, linestyle='--')
gl.xlabels_top = False ##关闭上侧坐标显示
gl.ylabels_right = False ##关闭右侧坐标显示
gl.xlabel_style={'size':5}
gl.ylabel_style={'size':5}
gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+0.1, 5))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+0.1, 5))

lynpatty

我自己顶一顶。。。

lynpatty

自己顶一顶

平凡之路

楼主不错，学习一下

meteorological

6666666666666666666