Satpy库绘制位置偏移

mushroom111

左侧是用下面代码画的，右图是micaps截图，对流位置偏北，咋解决呢（数据文件太大上传不了又咋办呢（累了

1. from satpy import Scene
   
2. from satpy import available_readers
   
3. from glob import glob
   
4. import matplotlib.pyplot as plt
   
5. import cartopy.crs as ccrs
   
6. import numpy as np
   
7. import os, glob
   
8. import datetime
   
9. from datetime import timedelta
   
10. import cartopy.io.shapereader as shpreader
    
11. # file = 'F:\\FY4A\\FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200517070000_20200517071459_0500M_V0001.HDF'
    
12. path = r'F:/file1/'
    
13. for index, filename in enumerate(os.listdir(path)):
    
14.     l1 = os.path.join(path, filename)
    
15.     filenames = glob.glob(l1)
    
16.     scn = Scene(filenames = filenames, reader = 'agri_fy4a_l1')
    
17.     scn.available_dataset_names()
    
18.     channel='C12'
    
19.     scn.load([channel])
    
20.     xmin,xmax,ymin,ymax=100,140,30,50
    
21.     scn2=scn.crop(ll_bbox=(xmin,ymin,xmax,ymax))
    
22.     # scn2.show(channel)
    
23.     fig=plt.figure(figsize=(12, 10), dpi = 600)
    
24.     area_def = scn2[channel].attrs['area']
    
25.     crs = area_def.to_cartopy_crs()
    
26.     datacrs=ccrs.PlateCarree()
    
27.     ax=plt.axes(projection = datacrs)# 添加地理特征
    
28.     # ax.coastlines(color='red', linewidth = 0.6)
    
29.     coastlines = "shp\\ne_50m_coastline.shp"  #预先准备好的50m海岸线地图
    
30.     adm2_shapes = list(shpreader.Reader(coastlines).geometries())    #读取地图文件信息
    
31.     ax.add_geometries(adm2_shapes, ccrs.PlateCarree(), edgecolor='red', linewidth = 0.6,  facecolor='none')
    
32.     ax.gridlines(crs=ccrs.PlateCarree(),draw_labels = {"bottom":"x","left":"y"},linestyle='--',alpha=0.5)# 显示图像
    
33.     img=ax.imshow(scn2[channel],extent = [xmin,xmax,ymin,ymax],
    
34.             origin='upper',cmap='Greys_r',transform=datacrs)
    
35.     cb=fig.colorbar(img,ax=ax,shrink=0.5,location="right",pad=0.02)
    
36.     ax.set_extent([xmin,xmax,ymin,ymax], crs=datacrs)
    
37.     cb.ax.tick_params(labelsize=8)
    
38.     # cb.set_label(label='Reflectance (%)', fontsize = 8)# 添加标题
    
39.     utc_time = filename[44:58]
    
40.     t_utc = datetime.datetime.strptime(utc_time, '%Y%m%d%H%M%S')
    
41.     t_bj = t_utc + timedelta(hours=8)
    
42.     ax.set_title('Channel12 ' + str(t_bj))
    
43. 
    
44.     # filename="F:\\FY4A_visC02_"+utc_time+".png"
    
45.     # plt.savefig(filename,dpi=600,bbox_inches='tight')
    
46.     # plt.close()

复制代码

edwardli

【气Py-47-卫星-哔哩哔哩】 https://b23.tv/TIrx4Rl