정보
-
업무명 : 파이썬 MOHID 해양순환모델 매핑 : 파일 읽기 및 함수 설정
-
작성자 : 박진만
-
작성일 : 2019-09-02
-
설 명 : MOHID 해양순환모델을 읽고 이를 매핑하기 위해 함수 및 basemap, 그리고 shp 파일을 로드하는 작업
-
수정이력 :
내용
[특징]
-
NC 파일을 읽고 이를 매핑하기 위해, shp, dbl 파일을 불러와 이를 구조화 한다.
[기능]
-
MOHID 해양순환모델 자료 읽기 (NetCDF 파일)
-
shp 파일 및 dbl 파일 읽기
-
기타 밑그림을 그리기 위한 작업 및 basemap 로드
[활용 자료]
-
MOHID 해양순환 모델, 부산지역 shp 파일 및 dbl 파일
[자료 처리 방안 및 활용 분석 기법]
-
없음
[사용법]
-
입력자료 (mohid nc 파일 및 shp, dbl 파일) 를 읽은 후 해당 소스코드를 실행한다.
[사용 OS]
-
Windows 10
[사용 언어]
-
Python v2.7
소스 코드
## Library 로드하기
import os
import pandas as pd
import numpy as np
import sys
import matplotlib.pyplot as plt
from dplython import *
# (DplyFrame, X, diamonds, select, sift, sample_n,
# sample_frac, head, arrange, mutate, group_by, summarize, DelayFunction)
from scipy.stats import linregress
from matplotlib import pyplot as plt
from IPython.display import Image
from mpl_toolkits.basemap import Basemap
from matplotlib.colors import Normalize
import matplotlib
import matplotlib.cm as cm
import seaborn as sns
from scipy.stats import linregress
from matplotlib import rcParams
from netCDF4 import Dataset
import struct
import binascii
from mpl_toolkits.basemap import addcyclic
from mpl_toolkits.basemap import shiftgrid
from netCDF4 import num2date, date2num, date2index
import datetime
from pyhdf.SD import SD, SDC
import h5py
import timeit
from datetime import datetime
from colormap import rgb2hex
from matplotlib.colors import LinearSegmentedColormap
from matplotlib import rc
from mpl_toolkits.axes_grid1 import make_axes_locatable
import matplotlib.patches as patches
from scipy.ndimage.filters import gaussian_filter
from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
from matplotlib.patches import PathPatch
import shapefile
#MATPLOTLIB 경고 무시
import warnings
warnings.filterwarnings("ignore")
#파일 인포 조회
path = "D:/ocean_model_out/mohid_hydr"
os.chdir(path)
os.getcwd()
ncfile_p = 'BS_MOHID_HYDR_2019021800.nc'
ncfile = Dataset(ncfile_p)
print(ncfile)
print(ncfile.variables.keys())
nc_value = ncfile.variables.keys()
#위,경도 추출
lat_staggered = ncfile.variables[u'lat_staggered'][:]
lon_staggered = ncfile.variables[ u'lon_staggered'][:]
#깊이, 마스크, ssh 추출
bathymetry = ncfile.variables[u'bathymetry'][:]
mask = ncfile.variables[u'mask'][:]
ssh = ncfile.variables[u'ssh'][:]
# U,V,VM 추출
U = ncfile.variables[u'u'][:]
V = ncfile.variables[u'v'][:]
VM = ncfile.variables[u'vm'][:]
# 염도, 수온 추출
salinity = ncfile.variables[u'salinity'][:]
temperature = ncfile.variables[u'temperature'][:]
# 마스크 포인트 추출
mask_OpenPoints = ncfile.variables[u'mask_OpenPoints'][:]
# 각 지역 POLYGON에 대한 파일 및 디렉토리 잡아주기
BS1_DIR = "D:/ocean_model_out/wave_ko/shpfile/BS_1Poly"
BS1_FILE = "BS_1Poly"
BS2_DIR = "D:/ocean_model_out/wave_ko/shpfile/BS_2Poly"
BS2_FILE = "BS_2Poly"
BS3_DIR = "D:/ocean_model_out/wave_ko/shpfile/BS_3Poly"
BS3_FILE = "BS_3Poly"
BS4_DIR = "D:/ocean_model_out/wave_ko/shpfile/BS_4Poly"
BS4_FILE = "BS_4Poly"
BSPOLY_DIR = "D:/ocean_model_out/wave_ko/shpfile/BSPoly"
BSPOLY_FILE = "BSPoly"
JPPOLY_DIR = "D:/ocean_model_out/wave_ko/shpfile/JPPoly"
JPPOLY_FILE = "JPPoly"
#;---------Domain Setting -------start- :: BS 1
# bs = 0 ;4 ;boundary space
# Ls = bs+110 ;bs+0 ;lat start
# Ms = bs+1 ;bs+0 ;lon start
# Lm = Ls+213 ;Lm-bs ;lat end
# Mm = Ms+240 ;Mm-bs ;lon end
#;--------------------------------end-
# lat = 660 lon = 836
################### BS #########################
minlon = min(lon)
maxlon = max(lon)
minlat = min(lat)
maxlat = max(lat)
################### BS1 ########################
bs1 = 0 #4 ;boundary space
Ls1 = bs1+110 #;bs+0 ;lat start
Ms1 = bs1+1 #;bs+0 ;lon start
Lm1 = Ls1+213 #;Lm-bs ;lat end
Mm1 = Ms1+240 #;Mm-bs ;lon end
################### BS2 ########################
bs2 = 0 #4 ;boundary space
Ls2 = bs2+210 #;bs+0 ;lat start
Ms2 = bs2+165 #;bs+0 ;lon start
Lm2 = Ls2+200 #;Lm-bs ;lat end
Mm2 = Ms2+240 #;Mm-bs ;lon end
################### BS3 ########################
bs3 = 0 #4 ;boundary space
Ls3 = bs3+335 #;bs+0 ;lat start
Ms3 = bs3+220 #;bs+0 ;lon start
Lm3 = Ls3+210 #;Lm-bs ;lat end
Mm3 = Ms3+240 #;Mm-bs ;lon end
################### BS4 ########################
bs4 = 0 #4 ;boundary space
Ls4 = bs4+440 #;bs+0 ;lat start
Ms4 = bs4+220 #;bs+0 ;lon start
Lm4 = Ls4+219 #;Lm-bs ;lat end
Mm4 = Ms4+300 #;Mm-bs ;lon end
print(Ls1,Lm1,Ms1,Mm1)
print(Ls2,Lm2,Ms2,Mm2)
print(Ls3,Lm3,Ms3,Mm3)
print(Ls4,Lm4,Ms4,Mm4)
MINLON_BS = min(lon)
MAXLON_BS = max(lon)
MINLAT_BS = min(lat)
MAXLAT_BS = max(lat)
MINLON_BS1 = lon[Ms1]
MAXLON_BS1 = lon[Mm1]
MINLAT_BS1 = lat[Ls1]
MAXLAT_BS1 = lat[Lm1]
MINLON_BS2 = lon[Ms2]
MAXLON_BS2 = lon[Mm2]
MINLAT_BS2 = lat[Ls2]
MAXLAT_BS2 = lat[Lm2]
MINLON_BS3 = lon[Ms3]
MAXLON_BS3 = lon[Mm3]
MINLAT_BS3 = lat[Ls3]
MAXLAT_BS3 = lat[Lm3]
MINLON_BS4 = lon[Ms4]
MAXLON_BS4 = lon[Mm4]
MINLAT_BS4 = lat[Ls4]
MAXLAT_BS4 = lat[Lm4]
bsmap = Basemap(projection='merc', lon_0=(MINLON_BS + MAXLON_BS) / 2.0 , lat_0=(MINLAT_BS + MAXLAT_BS) / 2.0,
llcrnrlon=MINLON_BS, llcrnrlat=MINLAT_BS,
urcrnrlon=MAXLON_BS, urcrnrlat=MAXLAT_BS,
resolution='c')
bsmap1 = Basemap(projection='merc', lon_0=(MINLON_BS1 + MAXLON_BS1) / 2.0 , lat_0=(MINLAT_BS1 + MAXLAT_BS1) / 2.0,
llcrnrlon=MINLON_BS1 , llcrnrlat=MINLAT_BS1,
urcrnrlon=MAXLON_BS1 , urcrnrlat=MAXLAT_BS1,
resolution='c')
bsmap2 = Basemap(projection='merc', lon_0=(MINLON_BS2 + MAXLON_BS2) / 2.0 , lat_0=(MINLAT_BS2 + MAXLAT_BS2) / 2.0,
llcrnrlon=MINLON_BS2, llcrnrlat=MINLAT_BS2,
urcrnrlon=MAXLON_BS2, urcrnrlat=MAXLAT_BS2,
resolution='c')
bsmap3 = Basemap(projection='merc', lon_0=(MINLON_BS3 + MAXLON_BS3) / 2.0 , lat_0=(MINLAT_BS3 + MAXLAT_BS3) / 2.0,
llcrnrlon=MINLON_BS3, llcrnrlat=MINLAT_BS3,
urcrnrlon=MAXLON_BS3, urcrnrlat=MAXLAT_BS3,
resolution='c')
bsmap4 = Basemap(projection='merc', lon_0=(MINLON_BS4 + MAXLON_BS4) / 2.0 , lat_0=(MINLAT_BS4 + MAXLAT_BS4) / 2.0,
llcrnrlon=MINLON_BS4, llcrnrlat=MINLAT_BS4,
urcrnrlon=MAXLON_BS4, urcrnrlat=MAXLAT_BS4,
resolution='c')
bsmap.readshapefile(BSPOLY_DIR, BSPOLY_FILE)
bsmap1.readshapefile(BS1_DIR, BS1_FILE)
bsmap2.readshapefile(BS2_DIR, BS2_FILE)
bsmap3.readshapefile(BS3_DIR, BS3_FILE)
bsmap4.readshapefile(BS4_DIR, BS4_FILE)
bsmap.readshapefile(JPPOLY_DIR, JPPOLY_FILE)
bsmap1.readshapefile(JPPOLY_DIR, JPPOLY_FILE)
bsmap2.readshapefile(JPPOLY_DIR, JPPOLY_FILE)
bsmap3.readshapefile(JPPOLY_DIR, JPPOLY_FILE)
bsmap4.readshapefile(JPPOLY_DIR, JPPOLY_FILE)
lon2D, lat2D = np.meshgrid(lon, lat)
lon2D_BS1, lat2D_BS1 = np.meshgrid(lon[Ms1:Mm1], lat[Ls1:Lm1])
lon2D_BS2, lat2D_BS2 = np.meshgrid(lon[Ms2:Mm2], lat[Ls2:Lm2])
lon2D_BS3, lat2D_BS3 = np.meshgrid(lon[Ms3:Mm3], lat[Ls3:Lm3])
lon2D_BS4, lat2D_BS4 = np.meshgrid(lon[Ms4:Mm4], lat[Ls4:Lm4])
BS_X, BS_Y = bsmap(lon2D,lat2D)
BS1_X, BS1_Y = bsmap1(lon2D_BS1, lat2D_BS1)
BS2_X, BS2_Y = bsmap2(lon2D_BS2, lat2D_BS2)
BS3_X, BS3_Y = bsmap3(lon2D_BS3, lat2D_BS3)
BS4_X, BS4_Y = bsmap4(lon2D_BS4, lat2D_BS4)
## bln 파일 로드
blnFilePath1 = "D:/ocean_model_out/mohid_hydr/blue-line.bln"
blnFilePath2 = "D:/ocean_model_out/mohid_hydr/deasan-habor-line.bln"
blnFilePath3 = "D:/ocean_model_out/mohid_hydr/green-line.bln"
blnFilePath4 = "D:/ocean_model_out/mohid_hydr/inchon-habor-line.bln"
blnFilePath5 = "D:/ocean_model_out/mohid_hydr/KMA-habor-line.bln"
blnFilePath6 = "D:/ocean_model_out/mohid_hydr/red-line.bln"
blue_line = pd.read_csv(blnFilePath1, sep='\s+', header=None)
deasan_habor_line = pd.read_csv(blnFilePath2, sep='\s+', header=None)
green_line = pd.read_csv(blnFilePath3, sep='\s+', header=None)
inchon_habor_line = pd.read_csv(blnFilePath4, sep='\s+', header=None)
KMA_habor_line = pd.read_csv(blnFilePath5, sep='\s+', header=None)
red_line = pd.read_csv(blnFilePath6, sep='\s+', header=None)
blue_linex = [[] for i in range(0,1000)]
blue_liney = [[] for i in range(0,1000)]
deasan_habor_linex = [[] for i in range(0,1000)]
deasan_habor_liney = [[] for i in range(0,1000)]
green_linex = [[] for i in range(0,1000)]
green_liney = [[] for i in range(0,1000)]
inchon_habor_linex = [[] for i in range(0,1000)]
inchon_habor_liney = [[] for i in range(0,1000)]
KMA_habor_linex = [[] for i in range(0,1000)]
KMA_habor_liney = [[] for i in range(0,1000)]
blue_linex = [[] for i in range(0,1000)]
blue_liney = [[] for i in range(0,1000)]
red_linex = [[] for i in range(0,1000)]
red_liney = [[] for i in range(0,1000)]
count_blue = 0
for i in range(0,len(blue_line)):
if( ((blue_line[0][i] - blue_line[1][i]) % 1) != 0 ):
blue_linex[count_blue].append(blue_line[0][i])
blue_liney[count_blue].append(blue_line[1][i])
else:
count_blue += 1;
count = 0
for i in range(0,len(deasan_habor_line)):
if( ((deasan_habor_line[0][i] - deasan_habor_line[1][i]) % 1) != 0 ):
deasan_habor_linex[count].append(deasan_habor_line[0][i])
deasan_habor_liney[count].append(deasan_habor_line[1][i])
else:
count += 1;
count_green = 0
for i in range(0,len(green_line)):
if( ((green_line[0][i] - green_line[1][i]) % 1) != 0 ):
green_linex[count_green].append(green_line[0][i])
green_liney[count_green].append(green_line[1][i])
else:
count_green += 1;
count = 0
for i in range(0,len(inchon_habor_line)):
if( ((inchon_habor_line[0][i] - inchon_habor_line[1][i]) % 1) != 0 ):
inchon_habor_linex[count].append(inchon_habor_line[0][i])
inchon_habor_liney[count].append(inchon_habor_line[1][i])
else:
count += 1;
count = 0
for i in range(0,len(KMA_habor_line)):
if( ((KMA_habor_line[0][i] - KMA_habor_line[1][i]) % 1) != 0 ):
KMA_habor_linex[count].append(KMA_habor_line[0][i])
KMA_habor_liney[count].append(KMA_habor_line[1][i])
else:
count += 1;
count_red = 0
for i in range(0,len(red_line)):
if( ((red_line[0][i] - red_line[1][i]) % 1) != 0 ):
red_linex[count_red].append(red_line[0][i])
red_liney[count_red].append(red_line[1][i])
else:
count_red += 1;
#[18,16,15,13,12]
#Depth 값에 따른 마스크 setting
bathymetry_0m_point = np.where(bathymetry>0,1.0,-999.0)
bathymetry_10m_point = np.where(bathymetry>10,1.0,-999.0)
bathymetry_25m_point = np.where(bathymetry>25,1.0,-999.0)
bathymetry_50m_point = np.where(bathymetry>50,1.0,-999.0)
bathymetry_75m_point = np.where(bathymetry>75,1.0,-999.0)
def switch(value,index,Ms=0,Mm=660,Ls=0,Lm=836):
if(index == 18):
value = value * bathymetry_0m_point[Ms:Mm,Ls:Lm]
elif(index == 16):
value = value * bathymetry_10m_point[Ms:Mm,Ls:Lm]
elif(index == 15):
value = value * bathymetry_25m_point[Ms:Mm,Ls:Lm]
elif(index == 13):
value = value * bathymetry_50m_point[Ms:Mm,Ls:Lm]
elif(index == 12):
value = value * bathymetry_75m_point[Ms:Mm,Ls:Lm]
else:
print("인덱스에 없는 value입니다. 처음의 value 값이 반환되었습니다.")
return(value)
# shp 파일 색칠하는 함수들..
def BSPOLY(map_object):
patches = []
for info, shape in zip(map_object.BSPoly_info, map_object.BSPoly):
patches.append( Polygon(np.array(shape), True) )
return(patches)
def BSPOLY1(map_object):
patches = []
for info, shape in zip(map_object.BS_1Poly_info, map_object.BS_1Poly):
patches.append( Polygon(np.array(shape), True) )
return(patches)
def BSPOLY2(map_object):
patches = []
for info, shape in zip(map_object.BS_2Poly_info, map_object.BS_2Poly):
patches.append( Polygon(np.array(shape), True) )
return(patches)
def BSPOLY3(map_object):
patches = []
for info, shape in zip(map_object.BS_3Poly_info, map_object.BS_3Poly):
patches.append( Polygon(np.array(shape), True) )
return(patches)
def BSPOLY4(map_object):
patches = []
for info, shape in zip(map_object.BS_4Poly_info, map_object.BS_4Poly):
patches.append( Polygon(np.array(shape), True) )
return(patches)
def JPPOLY(map_object):
patches = []
for info, shape in zip(map_object.JPPoly_info, map_object.JPPoly):
patches.append( Polygon(np.array(shape), True) )
return(patches)
def greenline(map_object):
for line in range(0,count_green,1):
#print(line)
aa, bb = map_object(green_linex[line][:],green_liney[line][:])
plt.plot(aa,bb, color='green')
return(0)
def redline(map_object):
for line in range(0,count_red,1):
#print(line)
aa, bb = map_object(red_linex[line][:],red_liney[line][:])
plt.plot(aa,bb, color='red')
return(0)
def blueline(map_object):
for line in range(0,count_blue,1):
#print(line)
aa, bb = map_object(blue_linex[line][:],blue_liney[line][:])
plt.plot(aa,bb, color='blue')
return(0)
# 시간 추출 (UTC and KST)
wet = 1104537600 #start time = 2005/01/01/00/00/00 UTC
print(lat.shape, lon.shape)
print(time)
timelen = len(time)
times = [0]*timelen
local_dt = [0]*timelen
print(times)
for i in range(0,timelen):
times[i] = datetime.utcfromtimestamp(time[int(i)] + wet).strftime('%Y-%m-%d %H') # to UTC
local_dt[i] = datetime.fromtimestamp(time[int(i)] + wet).strftime('%Y-%m-%d %H') # to KST
times = datetime.utcfromtimestamp(1558850400).strftime('%Y-%m-%d %H') # to UTC
# ncl color table
## 컬러바 만들기
# color_table = [[235,255,255,1],
# [235,255,255,1],
# [229,255,255],
# [211,255,255],
# [193,255,255],
# [175,243,255],
# [164,225,255],
# [157,207,255],
# [255,255,180],
# [255,255,160],
# [255,255,140],
# [255,255,120],
# [255,255,100],
# [255,255,50],
# [255,255,0],
# [255,202,108],
# [255,184,90],
# [255,166,72],
# [255,148,54],
# [255,130,36],
# [255,112,18],
# [255,110,0],
# [237,94,0],
# [255,72,72],
# [255,54,54],
# [255,18,18],
# [237,0,0],
# [219,0,0],
# [201,0,0],
# [165,0,0],
# [147,0,0],
# [129,0,0],
# [111,0,0],
# [93,0,0],
# [93,0,0] ];
color_table = [[235,255,255,1],
[229,255,255],
[211,255,255],
[193,255,255],
[175,243,255],
[164,225,255],
[157,207,255],
[255,255,180],
[255,255,160],
[255,255,140],
[255,255,120],
[255,255,100],
[255,255,50],
[255,255,0],
[255,202,108],
[255,184,90],
[255,166,72],
[255,148,54],
[255,130,36],
[255,112,18],
[255,110,0],
[237,94,0],
[255,72,72],
[255,54,54],
[255,18,18],
[237,0,0],
[219,0,0],
[201,0,0],
[165,0,0],
[147,0,0],
[129,0,0],
[111,0,0],
[93,0,0] ];
a = [0]*len(color_table)
for i in range(0, len(color_table)):
a[i] = rgb2hex(color_table[i][0],color_table[i][1],color_table[i][2])
b = ['#000000','#000000']
c = ["#00007F", "blue", "#007FFF", "cyan","#7FFF7F", "yellow", "#FF7F00", "red"]
cmap_jet = LinearSegmentedColormap.from_list('mycmap', a)
vcmap = LinearSegmentedColormap.from_list('vcmap', b)
cmap_jet1 = LinearSegmentedColormap.from_list('mycmap', c)
cmap_jet.set_under('w')
cmap_jet.set_over('#2F0000')
cmap_jet.set_bad('w')
cmap_jet1.set_under('#00007F')
cmap_jet1.set_over('#FF0000')
cmap_jet1.set_bad('w')
yy=np.arange(0,len(lat),20)
xx=np.arange(0,len(lon),20)
points=np.meshgrid(yy,xx)
cv_tick = np.linspace(0.0, 2.0, 31, endpoint=True)
cv_cbar_tick = np.linspace(0.0, 2.0, 11, endpoint=True)
tp_cbar_tick = np.linspace(0.0, 33.0, 12, endpoint=True)
tp_conf_tick = np.linspace(0.0, 33.0, 64, endpoint=True)
tp_cont_tick = np.linspace(0.0, 33.0, 34, endpoint=True)
sa_cbar_tick = np.linspace(25.0, 35.0, 6, endpoint=True)
sa_conf_tick = np.linspace(25.0, 35.0, 64, endpoint=True)
sa_cont_tick = np.linspace(25.0, 35.0, 21, endpoint=True)
참고문헌
[논문]
-
없음
[보고서]
- 없음
[URL]
- 없음
블로그에 대한 궁금하신 점을 문의하시면 자세히 답변드리겠습니다.
E. sangho.lee.1990@gmail.com & saimang0804@gmail.com
'프로그래밍 언어 > Python' 카테고리의 다른 글
| [Python] 파이썬 라디오미터 밝기온도 자료를 이용한 시계열 그래프 (0) | 2019.09.07 |
|---|---|
| [Python] 파이썬 과학 기술 계산 및 처리에 유용한 패키지 소개 (0) | 2019.09.05 |
| [Python] 파이썬 바이너리 (Binary) 형식인 천리안위성 1A호 (COMS/MI) 기상 위성 자료를 이용한 가시화 (0) | 2019.09.03 |
| [Python] 파이썬 아스키 (ASCII) 형식인 히마와리 8호 (Himawari-8/AHI) 및 Aqua/CERES 기상 위성 자료를 이용한 가시화 (0) | 2019.09.03 |
| [Python] 파이썬 MOHID 해양순환모델 결과 (남동해역 유속, 수온, 염분)를 이용한 가시화 (0) | 2019.09.03 |
최근댓글