작업26: 파이썬으로 대기 상단 연평균 일사량 계산하고 시각화하기1

정말 오랜만에 파이썬 글을 쓰네요..

 

시험 기간이기도 하고 연구실 구성하는거랑 이것 저것 일이 많아서 좀 걸렸네요 ㅠㅠ

 

그동안 했던 것들을 다시 하나씩 올려보려고 합니다

 

오늘은 대기 상단 연평균 일사량을 계산해보고 시각화해보려고 합니다

 

기후학이라는 과목에서 과제로 나왔습니다

 

다음에 기후학을 듣는 학생들도 도움이 됐으면 좋겠습니다

---

먼저 코드 공유전에 몇가지 공식들 설명이 필요합니다

포트란으로 적위와 이심률을 계산했던 포스팅이 있어서 보셨을만한 공식입니다

 

순서대로 태양의 적위, 감마, 이심률 공식입니다

 

포트란으로 계산해 둔 테이블이 있지만 이번에는 파이썬으로도 새로 계산해서 사용해보려고 합니다

 

다음은 시간각입니다

 

정확히는 일출 시간각을 구해서 사용합니다

 

이렇게 구해진 3가지 값을 이용해서 아래 공식을 계산합니다

대기 상단 일평균 일사량 공식입니다

 

여기서 파이는 위도입니다

 

위도의 범위는 -90º부터 90º까지 입니다

 

S는 태양상수로 1361입니다

 

바로 코드로 가보겠습니다

---

#일각 (gamma)
gam = (2 * np.pi * (np.arange(1, 366, 1) - 1))/365

#적위
de = []
for mn in gam:
  cor = (0.006918 - 0.399912*np.cos(mn)+ 0.0070257*np.sin(mn) - 0.006758*np.cos(2*mn) + 0.000907*np.sin(2*mn) - 0.002697*np.cos(3*mn) + 0.00148*np.sin(3*mn)) * (180/np.pi)
  de.append(cor)
  
#이심률
be = []
for r in np.arange(1, 366, 1):  
  Ec = 1 + 0.033 * np.cos((2 * np.pi * r)/365)
  be.append(Ec)

적위와 이심률은 계산되는 값의 갯수가 365개 입니다

 

각자 공식으로 계산을 해주고 빈 리스트에 넣어서 나중에 index를 부르는 방식으로 사용하기 위해

 

리스트로 만들어줬습니다

 

 

def cal_w(z):
  a = []

  dd1 = de

  for j in range(len(dd1)):
    tt = -np.tan(np.radians(z)) * np.tan(np.radians(dd1[j]))
    if tt > 1:
      w = 0
    elif tt < -1:
      w = np.radians(180)
    else:
      w = np.arccos(tt)

    a.append(w)
    zzz = ((np.array(a) * 180)/np.pi)
  globals()['a{}'.format(z)] = (zzz * 2 * np.pi)/360

다음은 시간각입니다

 

전체 코드는 함수 형태 입니다

 

위도별로 전부 계산을 해야하기 때문에 함수로 만들고 

 

나중에 np.range를 사용해서 한꺼번에 계산하는 방식으로 코드를 단순화시켰습니다

 

그래서 코드가 돌아가는 방식을 설명드리면 특정한 위도가 cal_w()안으로 들어가면

 

그 위도에서 1일 ~ 365일까지의 적위와 함께 365번 계산합니다

 

그렇게 때문에 위도별로 결과가 365개가 나와야합니다

 

만약 cal_w(60)이면 북위 60도에서의 1일 ~ 365일까지의 시간각이 총 365개가 나오는겁니다

---

시간각은 조건이 필요합니다

 

arccos을 사용하기 때문에 -tanπtanδ의 계산값이 arccos의 정의역을 넘어가면 계산이 되지 않습니다

 

그래서 만약에 -tanπtanδ가 -1보다 작은 값이 나오면 전부 π로 바꿔주고 1보다 크면 0으로 놓고 코드를 짜야합니다

 

  for j in range(len(dd1)):
    tt = -np.tan(np.radians(z)) * np.tan(np.radians(dd1[j]))
    if tt > 1:
      w = 0
    elif tt < -1:
      w = np.radians(180)
    else:
      w = np.arccos(tt)

그래서 코드가 이런 형태가 된겁니다

 

시간각 계산에는 적위가 사용되기 때문에 아까 말씀드린 것처럼 index를 불러오는 방식으로 dd1[j]를 사용했습니다

 

    a.append(w)
    zzz = ((np.array(a) * 180)/np.pi)
  globals()['a{}'.format(z)] = (zzz * 2 * np.pi)/360

그 후에 위에서 선언한 빈 리스트에 결과값을 하나씩 집어넣어줍니다

 

여기서 꽤 유용한 코드가 나옵니다

 

globals라는 코드입니다

 

지정해야할 변수 이름이 너무 많을 때 동적으로 변수를 만들어주는 코드입니다

 

a{}라는 형태로 이름을 만들 예정이고 format문을 활용해서 cal_w(60)의 결과가 a60이 될 수 있게 했습니다

 

이렇게 하면 변수를 일일이 지정하지 않아도 실행하는 만큼 생겨나게 됩니다

 

결과값은 각도가 아닌 라디안 형태로 바꿔주었습니다

 

π/180을 곱해도 되고 제가 한 것처럼 2π/360을 해주셔도 됩니다

 

저는 교수님이 2π/360으로 설명을 하셔서 저렇게 했습니다

for k in range(-90, 91, 1):
  cal_w(k)

실행은 이런식으로 위도별로 전부 실행해주었습니다

 

결과를 출력해서 확인해보면

잘 실행이 된 모습입니다

 

여기서 보시면 변수명 아래에 빨간줄이 있는 걸 보실 수 있는데

 

이는 globals로 생성한 변수들은 직접 선언한 변수들과는 다르게 전역 변수로

 

컴퓨터에 내적으로 저장시키는 것과 비슷하기 때문에 빨간줄이 생기는 겁니다

 

크게 신경 안쓰셔도 됩니다

---

이번 포스팅도 2편으로 나눠야할 거 같네요

 

글이 제 생각보다 많이 길어진 것 같습니다

 

다음 포스팅에서 이제 대기 상단 일사량 계산과 시각화를 마무리 해보겠습니다

 

긴 글 읽어주셔서 감사합니다!!