R語言、Python、MeteoInfo及Matlab繪製等值線圖

1、Matlab繪製等值線圖

註：南海小地圖請參考往期推文自行修改

clc;clear;close allf_hgt = 'F:/Rpython/lp8/ETOPO2v2g_f4.nc';lon=ncread(f_hgt,'x');lat=ncread(f_hgt,'y');h=ncread(f_hgt,'z');% plot出圖m_proj('Equidistant Cylindrical','long',[70 140],'lat',[15 55]);m_contourf(lon,lat,h',50,'linestyle','none');% 白化lon2 = [70 140];lat2 = [15 55];m_maskmap('F:/RMatlab/lp1/data/bou2_4p.shp',true,'lon',lon2,'lat',lat2,'m_map',true);hold onma=shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp'); m_line([ma(:).X],[ma(:).Y],'color','k');%繪製範圍內的地圖m_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[70:10:140],'ytick',[15:10:55]);c=colorbar('eastoutside','ticklength',0);caxis([0,8000])% 調用ncl色帶load('colorbar-mat/BkBlAqGrYeOrReViWh200_r.mat');colormap(BkBlAqGrYeOrReViWh200_r);% 調整colorbarax = gca;axpos = ax.Position;c.Position(3) = 0.5*c.Position(3);ax.Position = axpos;cbarrow;m_contour(lon,lat,h',[4000 6000],'ShowText','on','linestyle','-','linecolor','k');% 白化lon2 = [70 140];lat2 = [15 55];m_maskmap('F:/RMatlab/lp1/data/bou2_4p.shp',true,'lon',lon2,'lat',lat2,'m_map',true);hold onm_line([ma(:).X],[ma(:).Y],'color','k');%繪製範圍內的地圖m_grid('linestyle','none','box','fancy','fontsize',11,'tickdir','in','xtick',[70:10:140],'ytick',[15:10:55]);
2、R語言繪製等值線圖
註：其餘兩種白化方法請參考往期推文自行修改
library(ncdf4)library(ggplot2)library(ggspatial)bom <- nc_open("F:/Rpython/lp8/ETOPO2v2g_f4.nc")lon <- ncvar_get(bom, "x")lat <- ncvar_get(bom, "y")h <- ncvar_get(bom, "z")dimnames(h) <- list(lon, lat)library(tidyverse)library(RColorBrewer)library(reshape2)m <- h[,] %>%      melt(varnames = c("lon", "lat")) %>%      subset(!is.na(h))henan_shp <- "F:/Rpeng/32/data/henan.json"henan <- sf::read_sf(henan_shp)library(sf)m_df <-  st_as_sf(m,coords = c("lon", "lat"),crs = 4326)head(m_df)mask_region <- st_intersection(m_df,st_make_valid(henan))pal8 <- c("#33A02C", "#B2DF8A", "#FDBF6F", "#1F78B4", "#999999", "#E31A1C", "#E6E6E6", "#A6CEE3")henan <- ggplot()+geom_sf(data = mask_region,aes(color=value))+  geom_sf(data=henan,fill="NA",size=0.5,color="black")+theme_bw()+  scale_color_gradientn(colours=pal8,na.value = "transparent",  breaks = c(0,200,400,600,800,1000,1200,1400,1600,1800,2000),labels = c(0,200,400,600,800,1000,1200,1400,1600,1800,2000),limits = c(0,2000))+  theme(legend.key.height = grid::unit(2.6, "cm"))henan+annotation_scale(location = "bl") +  annotation_north_arrow(location="tr",which_north="false",style=north_arrow_fancy_orienteering)+  theme(axis.ticks.length = unit(0.2,"lines"), axis.ticks = element_line(color='black'),        axis.line = element_line(colour = "black"),axis.title.x=element_blank(),        axis.title.y=element_text(colour='black', size=16,face = "bold",vjust = 3),        axis.text.y=element_text(colour='black',size=16),        axis.text.x=element_text(colour = "black",size = 16,face = "bold",angle = 0,hjust = 0,vjust = 0.5),        plot.margin = margin(t = 5,r = 5,b = 5, l = 16, unit = "pt"),        text = element_text(colour = "black",size = 16,face = "bold"))+theme(panel.border = element_rect(fill=NA,color="black", size=1.1, linetype="solid"))+  geom_contour(data=m,aes(x=lon,y=lat,z=value),colour="white")
3、Python繪製等值線圖
注1：克裡金、反距離權重空間插值請參考往期推文自行修改
注2：南海小地圖請參考往期推文自行修改
import matplotlib.pyplot as pltimport numpy as npimport cartopy.crs as ccrsfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatterimport matplotlibimport gdalimport osfrom cartopy.io.shapereader import Readerfrom cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatterfrom matplotlib import rcParamsconfig={"font.family":'Times New Roman',"font.size":24,"mathtext.fontset":'stix'}rcParams.update(config)ds = gdal.Open('F:/Rpython/lp30/data/China_dem.tif')data = ds.ReadAsArray()gt = ds.GetGeoTransform()proj = ds.GetProjection()xres = gt[1]yres = gt[5]xmin = gt[0] + xres * 0.5xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5ymax = gt[3] - yres * 0.5ds = Nonexy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]proj= ccrs.PlateCarree()  fig=plt.figure(figsize=(16,12))extent=[70,140,15,55]  ax = fig.add_subplot(1, 1, 1,projection = proj)clevs=np.linspace(-200,7600,1000)map_car=ax.contourf(xy_source[0,::],xy_source[1,::],data.T,np.arange(-200,7600,1000),transform=ccrs.PlateCarree(),cmap='gist_rainbow',extend='both')map_car2=ax.contour(xy_source[0,::],xy_source[1,::],data.T,levels=[3500],colors='k')b=plt.colorbar(map_car,shrink=0.66,orientation='vertical',pad=0.035,aspect=20,extend='both')ax.set_extent(extent,crs = ccrs.PlateCarree())ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))ax.yaxis.set_major_formatter(LatitudeFormatter())ax.set_xticks(np.arange(70,140.1,10),crs=ccrs.PlateCarree())ax.set_yticks(np.arange(15,55.1,10),crs=ccrs.PlateCarree())SHP = 'F:/Rpython/lp27/data/'ax.add_geometries(Reader(os.path.join(SHP,'china2.shp')).geometries(),ccrs.PlateCarree(),facecolor='none',edgecolor='k',linewidth=1)ax0 = plt.gca()   ax0.outline_patch.set_linewidth(1)    plt.savefig('F:/Rpython/lp30/plot26.png',dpi=600)plt.show()
4、MeteoInfo繪製等值線圖
註：其餘參數設定請參考官網自行修改
fn = 'F:/RMeteoInfo/data/data2/dixing-china.nc'f = addfile(fn)data = f['z'][:,:] ax=axesm(tickfontsize=16,axis=True)world = shaperead('F:/RMeteoInfo/data/map/country1.shp') lchina = shaperead('F:/RMeteoInfo/data/map/bou2_4p.shp')lchina2 = shaperead('F:/RMeteoInfo/data/map/bou2_4l.shp')geoshow(lchina,edgecolor='k',size=0.8)bou1_layer=geoshow(lchina2,edgecolor='k',size=0.6)levels = arange(0, 8000, 1000)layer = imshowm(data,levels,cmap='temp_diff_18lev')masklayer(lchina, [layer])layer2 = contour(data,color='k',smooth=False)masklayer(lchina, [layer2])colorbar(layer,orientation='horizontal',ticklen=0,extendrect=False,shrink=1,aspect=50,label='Height(m)',bold=True,fontsize=15)yaxis(tickvisible=True,location='left',tickwidth=2,linewidth=2,ticklength=3) yaxis(tickvisible=False,location='right',tickwidth=2,linewidth=2,ticklength=4) xaxis(tickvisible=False,location='top',tickwidth=2,linewidth=2,ticklength=4) xaxis(tickvisible=True,location='bottom',tickwidth=2,linewidth=2,ticklength=3)xlim(70,140)ylim(15,55)xticks(arange(70, 142, 10),bold=True,fontsize=15)yticks(arange(15, 56, 10),bold=True,fontsize=15)xlabel('Longitude',bold=True,fontsize=15)ylabel('Latitude',bold=True,fontsize=15)ax.scale_bar(0.3,0.38,width=200,linewidth=1,bold=True,fontsize=12,bartype='alternating_bar')ax.north_arrow(0.3,0.9,width=40,height=40,linewidth=1,bold=True,fontsize=12)axesm(position=[0.648,0.273,0.15,0.2], axison=False, frameon=True)geoshow(lchina,edgecolor='k',size=0.8)bou1_layer=geoshow(lchina2,edgecolor='k',size=0.6)levels = arange(0, 8000, 1000)layer = imshowm(data,levels,cmap='temp_diff_18lev')masklayer(lchina, [layer])xlim(106, 123)ylim(2, 23)savefig('F:/RMeteoInfo/test/plot90.png',dpi=800)