Skip to content

Instantly share code, notes, and snippets.

@MeteoBoy4 MeteoBoy4/nc2wrfda.py
Last active Oct 5, 2017

Embed
What would you like to do?
Preparing precipitation data for WRFDA 4DVAR
#! /home/meteoboy4/anaconda/bin/python
import write_rainobs
import os
import numpy as np
import xarray as xr
import pandas as pd
MISSING_R = -888888.
dirvar = '/run/media/MeteoBoy4/Data/MData/TRMM'
data = 'TP.nc'
ncfile = xr.open_dataset(os.path.join(dirvar, data))
latitude = ncfile['LAT']
longitude = ncfile['LON']
times = ncfile['TIME']
rain_all = ncfile['PCP'].values
lon_matrix, lat_matrix = np.meshgrid(longitude, latitude)
lon = lon_matrix.flatten(order='C')
lat = lat_matrix.flatten(order='C')
N = lon.size
assert len(times) % 2 == 1
for i in range(1, len(times), 2):
rain_one = 6 * (rain_all[i] + rain_all[i+1]) / 2.
rain = rain_one.flatten(order='C')
rain[np.isnan(rain)] = MISSING_R
assert np.all(np.logical_or(rain>=0.0, rain<=MISSING_R))
date = pd.to_datetime(str(times[i+1].values)).strftime('%Y%m%d%H')
suffix = date + '.06h'
date_time = pd.to_datetime(str(times[i+1].values-1)).strftime('%Y-%m-%d_%H:%M:%S')
print(suffix)
print(date_time)
if i == 11:
write_rainobs.writerainobs(lat, lon, date_time, MISSING_R, suffix, rain)
!-*- f90 -*- - file contains F90 code in free format
SUBROUTINE WRITERAINOBS (N, LAT, LON, DATE, MISSING_R, SUFFIX, RAIN)
IMPLICIT NONE
!f2py intent(in) :: N, MISSING_R
!f2py intent(in) :: LAT, LON
!f2py intent(in) :: DATE, SUFFIX
!f2py intent(in) :: RAIN
INTEGER, INTENT(IN) :: N
REAL, INTENT(IN) :: LAT(N)
REAL, INTENT(IN) :: LON(N)
REAL, INTENT(IN) :: MISSING_R
CHARACTER(LEN=19), INTENT(IN) :: DATE
CHARACTER(LEN=14), INTENT(IN) :: SUFFIX
REAL, INTENT(IN) :: RAIN(N)
INTEGER :: I
CHARACTER(LEN=80) :: FILENAME
CHARACTER(LEN=120) :: FMT_INFO, &
FMT_EACH
LOGICAL :: CONNECTED
! 1. FORMAT
! ---------
!INFO = PLATFORM, DATE, LEVELS, LATITUDE, LONGITUDE, ELEVATION, ID.
!EACH = HEIGHT, RAINFALL DATA, QC, ERROR
FMT_INFO = '(A12, 1X, A19, 1X, I6, 2(F12.3,2X), F8.1, 1X, A5)'
FMT_EACH = '(F12.3, F12.3, I4, F12.3)'
! 2. OPEN A FILE FOR OBSERVATIONS OUTPUT
! --------------------------------------
! 2.1 OUTPUT FILE
! ---------------
WRITE(FILENAME, FMT = '(A, A)') 'ob.rain.', SUFFIX
WRITE(0, '(A)') 'Write precipitation data '
INQUIRE (UNIT = 999, OPENED = CONNECTED)
IF (.NOT. CONNECTED) THEN
OPEN (UNIT= 999, &
FILE = FILENAME, &
FORM = 'FORMATTED', &
ACCESS = 'SEQUENTIAL', &
STATUS = 'REPLACE')
END IF
REWIND (UNIT = 999)
! 3. FILE HEADER AND DATA
! -----------------------
! 3.1 TOTAL NUMBER OF OBSERVATIONS CONTAINED IN FILE
! ---------------------------------------------------
WRITE(UNIT = 999, FMT = '((A,I7),A)', ADVANCE = 'no') &
"TOTAL =", N,","
! 3.2 MISSING VALUE FLAG
! ------------------------
WRITE(UNIT = 999, FMT = '((A, F8.0), A)') "MISS. =", MISSING_R,","
! 3.3 FORMAT
! ------------
WRITE(UNIT = 999, FMT = '(A)') &
"INFO = PLATFORM, DATE, LEVELS, LATITUDE, LONGITUDE, ELEVATION, ID."
WRITE(UNIT = 999, FMT = '(A)' ) &
"EACH = HEIGHT, RAINFALL DATA, QC, ERROR"
! 3.4 DATA
! --------
WRITE (UNIT = 999, FMT = '(A)') &
"#------------------------------------------------------------------------------#"
DO I=1,N
WRITE(UNIT = 999, FMT = TRIM(FMT_INFO)) &
"FM-129 RAIN ", DATE, 1, LAT(I), LON(I),99.9, "12345"
WRITE(UNIT = 999, FMT = TRIM(FMT_EACH)) &
99.9, RAIN(I), 88, 2.000
END DO
CLOSE(999)
END SUBROUTINE WRITERAINOBS
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.