Python geojson文件可视化

最新推荐文章于 2025-03-07 14:49:43 发布

十光@夏天

最新推荐文章于 2025-03-07 14:49:43 发布

阅读量1.1k

点赞数 1

文章标签： python

本文链接：https://blog.csdn.net/weixin_45577825/article/details/120446017

版权

本文通过Python编程实现气象数据的读取与处理，并采用地理信息系统技术完成边界层高度变化的可视化展示，对比了不同日期间的变化情况。

摘要由CSDN通过智能技术生成

# 从整体上来看，先看边界层高度的情况
# 整理气象数据
# 读数据
import netCDF4 as nc
import pandas as pd
import numpy as np

file1 = r'.\adaptor.mars.internal-1631981155.7723448-28568-14-37c609d7-a545-4a66-b4e3-3f73bf08819b.nc'
dataset1 = nc.Dataset(file1)
lon1=dataset1.variables['longitude'][:]
lat1=dataset1.variables['latitude'][:]
BLH=dataset1.variables['blh'][:]

BLH20190203=np.mean(BLH[745:2160,:,:],axis=0)
BLH20200203=np.mean(BLH[3625:5064,:,:],axis=0)
delta=BLH20190203-BLH20200203

import geopandas
from shapely import geometry
import matplotlib.pyplot as plt
import json

data = geopandas.read_file('中国.json')
a=data.geometry.to_json()
polygon_dict = json.loads(a)

features=len(polygon_dict["features"][:])
point_list=[]
for i in range(features):
    #print('features****************************************')
    coordinates=len(polygon_dict["features"][i]["geometry"]["coordinates"][:])
    for j in range(coordinates):
        #print('coordinates-----------------------------------')
        rows=int(len(np.array(polygon_dict["features"][i]["geometry"]["coordinates"][j]).flatten())/2)
        
        point=(np.array(polygon_dict["features"][i]["geometry"]["coordinates"][j]).flatten()).reshape(rows,2)
        #print(point)
        point_list.append(point)
ini=point_list[0]
num=len(point_list)
for k in range(1,num):
    ini=np.vstack((ini,point_list[k]))

import matplotlib.pyplot as plt
import matplotlib.path as mpath
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.mpl.ticker as cticker
import cartopy.io.shapereader as shpreader
import seaborn as sns
from matplotlib import rcParams
import matplotlib as mpl
mpl.rcParams["font.family"] = 'Times New Roman'  #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 16  #字体大小

font = {'family' : 'Times New Roman',
        'color'  : 'black',
        'weight' : 'normal',
        'size'   : 25,
        }
fig1 = plt.figure(figsize=(14,7))

########################################################################################另外三幅子图
ax2 = fig1.add_subplot(1,1,1, projection=ccrs.PlateCarree())
gl=ax2.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                  linewidth=1, color='gray', alpha=0.5, linestyle='--')

ax2.add_feature(cfeature.COASTLINE.with_scale('110m'))
img_extent1=[107,117,28,34]
ax2.set_extent(img_extent1, ccrs.PlateCarree())

lon,lat=np.meshgrid(lon1,lat1)
c7 = ax2.pcolormesh(lon,lat,-delta,transform=ccrs.PlateCarree(),cmap=plt.cm.bwr,vmin=-70,vmax=70,alpha=1)
cb=plt.colorbar(c7)
cb.ax.tick_params(labelsize=16)  #设置色标刻度字体大小
ax2.scatter(ini[:,0],ini[:,1],transform=ccrs.PlateCarree(),marker='.',color='k')

plt.savefig('test.jpg')

在这里插入图片描述