forked from 3rdparty/wrf-python
Browse Source
Also reorganized the test directory. Various snippets that may or may not be of use anymore added the the misc directory. Misc. NCL scripts used during development added to ncl directory.lon0
Bill Ladwig
6 years ago
31 changed files with 1709 additions and 580 deletions
@ -0,0 +1,40 @@ |
|||||||
|
|
||||||
|
from netCDF4 import Dataset |
||||||
|
|
||||||
|
VARS_TO_KEEP = ("Times", "XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", |
||||||
|
"XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", |
||||||
|
"T2", "PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD", |
||||||
|
"QGRAUP", "QRAIN", "QSNOW", "QICE", "MAPFAC_M", "MAPFAC_U", |
||||||
|
"MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC", |
||||||
|
"I_RAINC", "I_RAINNC") |
||||||
|
|
||||||
|
INPUT_FILE = "wrfout_d02_2005-08-28_00:00:00" |
||||||
|
OUTPUT_FILE = "wrfout_problem_file" |
||||||
|
DESIRED_TIME_INDEX = 0 |
||||||
|
|
||||||
|
with Dataset(INPUT_FILE) as infile, Dataset(OUTPUT_FILE, mode="w") as outfile: |
||||||
|
# Copy the global attributes |
||||||
|
outfile.setncatts(infile.__dict__) |
||||||
|
|
||||||
|
# Copy Dimensions |
||||||
|
for name, dimension in infile.dimensions.items(): |
||||||
|
if name != "Time": |
||||||
|
dimsize = len(dimension) |
||||||
|
outfile.createDimension(name, dimsize) |
||||||
|
else: |
||||||
|
outfile.createDimension(name, 1) |
||||||
|
|
||||||
|
# Copy Variables |
||||||
|
for name, variable in infile.variables.iteritems(): |
||||||
|
if name not in VARS_TO_KEEP: |
||||||
|
continue |
||||||
|
|
||||||
|
outvar = outfile.createVariable(name, variable.datatype, |
||||||
|
variable.dimensions, zlib=True) |
||||||
|
|
||||||
|
if len(variable.dimensions) > 1: |
||||||
|
outvar[:] = variable[DESIRED_TIME_INDEX, :] |
||||||
|
else: |
||||||
|
outvar[:] = variable[DESIRED_TIME_INDEX] |
||||||
|
|
||||||
|
outvar.setncatts(variable.__dict__) |
||||||
@ -0,0 +1,70 @@ |
|||||||
|
from __future__ import print_function, division |
||||||
|
|
||||||
|
import os |
||||||
|
import numpy as np |
||||||
|
from netCDF4 import Dataset |
||||||
|
from wrf import getvar, ALL_TIMES, to_np |
||||||
|
import xarray |
||||||
|
|
||||||
|
filename_list = ["/Users/ladwig/Documents/wrf_files/" |
||||||
|
"wrf_vortex_single/wrfout_d02_2005-08-28_00:00:00", |
||||||
|
"/Users/ladwig/Documents/wrf_files/wrf_vortex_single/" |
||||||
|
"wrfout_d02_2005-08-28_03:00:00", |
||||||
|
"/Users/ladwig/Documents/wrf_files/wrf_vortex_single/" |
||||||
|
"wrfout_d02_2005-08-28_06:00:00", |
||||||
|
"/Users/ladwig/Documents/wrf_files/wrf_vortex_single/" |
||||||
|
"wrfout_d02_2005-08-28_09:00:00"] |
||||||
|
|
||||||
|
result_shape = (4, 1, 96, 96) |
||||||
|
|
||||||
|
# Let's get the first time so we can copy the metadata later |
||||||
|
f = Dataset(filename_list[0]) |
||||||
|
# By setting squeeze to False, you'll get all the dimension names. |
||||||
|
z1 = getvar(f, "T2", 0, squeeze=False) |
||||||
|
xlat = getvar(f, "XLAT", 0) |
||||||
|
xlong = getvar(f, "XLONG", 0) |
||||||
|
|
||||||
|
|
||||||
|
z_final = np.empty(result_shape, np.float32) |
||||||
|
|
||||||
|
# Modify this number if using more than 1 time per file |
||||||
|
times_per_file = 1 |
||||||
|
|
||||||
|
data_times = [] |
||||||
|
xtimes = [] |
||||||
|
for timeidx in range(result_shape[0]): |
||||||
|
# Compute the file index and the time index inside the file |
||||||
|
fileidx = timeidx // times_per_file |
||||||
|
file_timeidx = timeidx % times_per_file |
||||||
|
|
||||||
|
f = Dataset(filename_list[fileidx]) |
||||||
|
z = getvar(f, "T2", file_timeidx) |
||||||
|
t = getvar(f, "Times", file_timeidx) |
||||||
|
xt = getvar(f, "xtimes", file_timeidx) |
||||||
|
data_times.append(to_np(t)) |
||||||
|
xtimes.append(to_np(xt)) |
||||||
|
z_final[timeidx, :] = z[:] |
||||||
|
f.close() |
||||||
|
|
||||||
|
# Let's make the metadata. Dimension names should copy easily if you set |
||||||
|
# sqeeze to False, otherwise you can just set them yourself is a tuple of |
||||||
|
# dimension names. Since you wanted |
||||||
|
# to keep the bottom_top dimension for this 2D variable (which is normally |
||||||
|
# removed), I'm doing this manually. |
||||||
|
z_dims = ["Time", "bottom_top", "south_north", "west_east"] |
||||||
|
|
||||||
|
# Xarray doesn't copy coordinates easily (it always complains about shape |
||||||
|
# mismatches), so do this manually |
||||||
|
z_coords = {} |
||||||
|
z_coords["Time"] = data_times |
||||||
|
z_coords["XTIME"] = ("Time",), xtimes |
||||||
|
z_coords["XLAT"] = ("south_north", "west_east"), xlat |
||||||
|
z_coords["XLONG"] = ("south_north", "west_east"), xlong |
||||||
|
z_name = "T2" |
||||||
|
|
||||||
|
# Attributes copy nicely |
||||||
|
z_attrs = {} |
||||||
|
z_attrs.update(z1.attrs) |
||||||
|
|
||||||
|
z_with_meta = xarray.DataArray(z_final, coords=z_coords, dims=z_dims, |
||||||
|
attrs=z_attrs, name=z_name) |
||||||
@ -0,0 +1,58 @@ |
|||||||
|
from netCDF4 import Dataset |
||||||
|
import numpy as np |
||||||
|
import matplotlib.pyplot as plt |
||||||
|
from matplotlib.cm import get_cmap |
||||||
|
import cartopy.crs as crs |
||||||
|
from cartopy.feature import NaturalEarthFeature |
||||||
|
|
||||||
|
from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy, |
||||||
|
cartopy_xlim, cartopy_ylim) |
||||||
|
|
||||||
|
# Open the NetCDF file |
||||||
|
ncfile = Dataset("wrfout_d01_2016-10-07_00_00_00") |
||||||
|
|
||||||
|
# Extract the pressure and wind variables |
||||||
|
|
||||||
|
p = getvar(ncfile, "pressure") |
||||||
|
# Note: This is m/s. |
||||||
|
ua = getvar(ncfile, "ua") |
||||||
|
va = getvar(ncfile, "va") |
||||||
|
|
||||||
|
# Interpolate u, and v winds to 950 hPa |
||||||
|
u_950 = interplevel(ua, p, 950) |
||||||
|
v_950 = interplevel(va, p, 950) |
||||||
|
|
||||||
|
# Get the lat/lon coordinates |
||||||
|
lats, lons = latlon_coords(u_950) |
||||||
|
|
||||||
|
# Get the map projection information |
||||||
|
cart_proj = get_cartopy(u_950) |
||||||
|
|
||||||
|
# Create the figure |
||||||
|
fig = plt.figure(figsize=(12, 9)) |
||||||
|
ax = plt.axes(projection=cart_proj) |
||||||
|
|
||||||
|
# Download and add the states and coastlines |
||||||
|
states = NaturalEarthFeature(category='cultural', scale='50m', |
||||||
|
facecolor='none', |
||||||
|
name='admin_1_states_provinces_shp') |
||||||
|
ax.add_feature(states, linewidth=0.5) |
||||||
|
ax.coastlines('50m', linewidth=0.8) |
||||||
|
|
||||||
|
# Add the 950 hPa wind barbs, only plotting every 'thin'ed barb. Adjust thin |
||||||
|
# as needed. Also, no scaling is done for the arrows, so you might need to |
||||||
|
# mess with the scale argument. |
||||||
|
thin = 50 |
||||||
|
plt.quiver(to_np(lons[::thin, ::thin]), to_np(lats[::thin, ::thin]), |
||||||
|
to_np(u_950[::thin, ::thin]), to_np(v_950[::thin, ::thin]), |
||||||
|
transform=crs.PlateCarree()) |
||||||
|
|
||||||
|
# Set the map bounds |
||||||
|
ax.set_xlim(cartopy_xlim(u_950)) |
||||||
|
ax.set_ylim(cartopy_ylim(v_950)) |
||||||
|
|
||||||
|
ax.gridlines() |
||||||
|
|
||||||
|
plt.title("Arrows") |
||||||
|
|
||||||
|
plt.show() |
||||||
@ -1,30 +1,29 @@ |
|||||||
|
|
||||||
import matplotlib.pyplot as plt |
import matplotlib.pyplot as plt |
||||||
from mpl_toolkits.basemap import Basemap |
from mpl_toolkits.basemap import Basemap |
||||||
|
|
||||||
|
|
||||||
def main(): |
def main(): |
||||||
bm = Basemap(projection = "rotpole", |
bm = Basemap(projection="rotpole", |
||||||
o_lat_p = 36.0, |
o_lat_p=36.0, |
||||||
o_lon_p = 180.0, |
o_lon_p=180.0, |
||||||
llcrnrlat = -10.590603, |
llcrnrlat=-10.590603, |
||||||
urcrnrlat = 46.591976, |
urcrnrlat=46.591976, |
||||||
llcrnrlon = -139.08585, |
llcrnrlon=-139.08585, |
||||||
urcrnrlon = 22.661009, |
urcrnrlon=22.661009, |
||||||
lon_0 = -106.0, |
lon_0=-106.0, |
||||||
rsphere = 6370000, |
rsphere=6370000, |
||||||
resolution = 'l') |
resolution='l') |
||||||
|
|
||||||
fig = plt.figure(figsize=(8,8)) |
fig = plt.figure(figsize=(8, 8)) |
||||||
ax = fig.add_axes([0.1,0.1,0.8,0.8]) |
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) |
||||||
|
|
||||||
bm.drawcoastlines(linewidth=.5) |
bm.drawcoastlines(linewidth=.5) |
||||||
|
|
||||||
print bm.proj4string |
print(bm.proj4string) |
||||||
|
|
||||||
plt.savefig("basemap_map.png") |
plt.savefig("basemap_map.png") |
||||||
plt.close(fig) |
plt.close(fig) |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
if __name__ == "__main__": |
if __name__ == "__main__": |
||||||
main() |
main() |
||||||
|
|||||||
@ -1,26 +1,25 @@ |
|||||||
|
# Test snippet for f2py |
||||||
import numpy as np |
import numpy as np |
||||||
|
|
||||||
import wrf._wrffortran |
import wrf._wrffortran |
||||||
errlen = int(wrf._wrffortran.constants.errlen) |
errlen = int(wrf._wrffortran.constants.errlen) |
||||||
|
|
||||||
|
a = np.ones((3, 3, 3)) |
||||||
a = np.ones((3,3,3)) |
b = np.zeros((3, 3, 3, 3)) |
||||||
b = np.zeros((3,3,3,3)) |
|
||||||
c = np.zeros(errlen, "c") |
c = np.zeros(errlen, "c") |
||||||
errstat = np.array(0) |
errstat = np.array(0) |
||||||
errmsg = np.zeros(errlen, "c") |
errmsg = np.zeros(errlen, "c") |
||||||
|
|
||||||
c[:] = "Test String" |
c[:] = "Test String" |
||||||
|
|
||||||
|
|
||||||
for i in xrange(2): |
for i in xrange(2): |
||||||
outview = b[i,:] |
outview = b[i, :] |
||||||
outview = outview.T |
outview = outview.T |
||||||
q = wrf._wrffortran.testfunc(a,outview,c,errstat=errstat,errmsg=errmsg) |
q = wrf._wrffortran.testfunc(a, outview, c, errstat=errstat, |
||||||
print errstat |
errmsg=errmsg) |
||||||
|
print(errstat) |
||||||
|
|
||||||
|
|
||||||
str_bytes = (bytes(char).decode("utf-8") for char in errmsg[:]) |
str_bytes = (bytes(char).decode("utf-8") for char in errmsg[:]) |
||||||
print repr(errmsg) |
print(repr(errmsg)) |
||||||
print "".join(str_bytes).strip() |
print("".join(str_bytes).strip()) |
||||||
@ -0,0 +1,230 @@ |
|||||||
|
# Hastily made script to read WPS intermediate files |
||||||
|
from __future__ import print_function |
||||||
|
import struct |
||||||
|
|
||||||
|
import numpy as np |
||||||
|
|
||||||
|
# Number of bytes used at the start and end of a fortran record to |
||||||
|
# indicate the record size |
||||||
|
SIZE_BYTES = 4 |
||||||
|
|
||||||
|
|
||||||
|
class WPSData(object): |
||||||
|
def __init__(self, ifv=None, hdate=None, xfcst=None, map_source=None, |
||||||
|
field=None, units=None, desc=None, xlvl=None, nx=None, |
||||||
|
ny=None, iproj=None, startloc=None, startlat=None, |
||||||
|
startlon=None, deltalat=None, deltalon=None, nlats=None, |
||||||
|
dx=None, dy=None, xlonc=None, truelat1=None, truelat2=None, |
||||||
|
earth_radius=None, is_wind_earth_rel=None, slab=None): |
||||||
|
|
||||||
|
self.ifv = ifv |
||||||
|
self.hdate = hdate |
||||||
|
self.xfcst = xfcst |
||||||
|
self.map_source = map_source |
||||||
|
self.field = field |
||||||
|
self.units = units |
||||||
|
self.desc = desc |
||||||
|
self.xlvl = xlvl |
||||||
|
self.nx = nx |
||||||
|
self.ny = ny |
||||||
|
self.iproj = iproj |
||||||
|
self.startloc = startloc |
||||||
|
self.startlat = startlat |
||||||
|
self.startlon = startlon |
||||||
|
self.deltalat = deltalat |
||||||
|
self.deltalon = deltalon |
||||||
|
self.nlats = nlats |
||||||
|
self.dx = dx |
||||||
|
self.dy = dy |
||||||
|
self.xlonc = xlonc |
||||||
|
self.truelat1 = truelat1 |
||||||
|
self.truelat2 = truelat2 |
||||||
|
self.earth_radius = earth_radius |
||||||
|
self.is_wind_earth_rel = is_wind_earth_rel |
||||||
|
self.slab = slab |
||||||
|
|
||||||
|
|
||||||
|
def _parse_record1(data): |
||||||
|
result = {} |
||||||
|
result["ifv"] = struct.unpack(">i", data) |
||||||
|
|
||||||
|
return result |
||||||
|
|
||||||
|
|
||||||
|
def _parse_record2(data): |
||||||
|
result = {} |
||||||
|
parsed = struct.unpack(">24sf32s9s25s46sfiii", data) |
||||||
|
result["hdate"] = parsed[0] |
||||||
|
result["xfcst"] = parsed[1] |
||||||
|
result["map_source"] = parsed[2] |
||||||
|
result["field"] = parsed[3] |
||||||
|
result["units"] = parsed[4] |
||||||
|
result["desc"] = parsed[5] |
||||||
|
result["xlvl"] = parsed[6] |
||||||
|
result["nx"] = parsed[7] |
||||||
|
result["ny"] = parsed[8] |
||||||
|
result["iproj"] = parsed[9] |
||||||
|
|
||||||
|
return result |
||||||
|
|
||||||
|
|
||||||
|
def _parse_record3(data, iproj): |
||||||
|
result = {} |
||||||
|
if iproj == 0: |
||||||
|
fmt = ">8sfffff" |
||||||
|
parsed = struct.unpack(fmt, data) |
||||||
|
result["startloc"] = parsed[0] |
||||||
|
result["startlat"] = parsed[1] |
||||||
|
result["startlon"] = parsed[2] |
||||||
|
result["deltalat"] = parsed[3] |
||||||
|
result["deltalon"] = parsed[4] |
||||||
|
result["earth_radius"] = parsed[5] |
||||||
|
elif iproj == 1: |
||||||
|
fmt = ">8sffffff" |
||||||
|
parsed = struct.unpack(fmt, data) |
||||||
|
result["startloc"] = parsed[0] |
||||||
|
result["startlat"] = parsed[1] |
||||||
|
result["startlon"] = parsed[2] |
||||||
|
result["dx"] = parsed[3] |
||||||
|
result["dy"] = parsed[4] |
||||||
|
result["truelat1"] = parsed[5] |
||||||
|
result["earth_radius"] = parsed[6] |
||||||
|
elif iproj == 3: |
||||||
|
fmt = ">8sffffffff" |
||||||
|
parsed = struct.unpack(fmt, data) |
||||||
|
result["startloc"] = parsed[0] |
||||||
|
result["startlat"] = parsed[1] |
||||||
|
result["startlon"] = parsed[2] |
||||||
|
result["dx"] = parsed[3] |
||||||
|
result["dy"] = parsed[4] |
||||||
|
result["xlonc"] = parsed[5] |
||||||
|
result["truelat1"] = parsed[6] |
||||||
|
result["truelat2"] = parsed[7] |
||||||
|
result["earth_radius"] = parsed[8] |
||||||
|
elif iproj == 4: |
||||||
|
fmt = ">8sfffff" |
||||||
|
parsed = struct.unpack(fmt, data) |
||||||
|
result["startloc"] = parsed[0] |
||||||
|
result["startlat"] = parsed[1] |
||||||
|
result["startlon"] = parsed[2] |
||||||
|
result["nlats"] = parsed[3] |
||||||
|
result["deltalon"] = parsed[4] |
||||||
|
result["earth_radius"] = parsed[5] |
||||||
|
elif iproj == 5: |
||||||
|
fmt = ">8sfffffff" |
||||||
|
parsed = struct.unpack(fmt, data) |
||||||
|
result["startloc"] = parsed[0] |
||||||
|
result["startlat"] = parsed[1] |
||||||
|
result["startlon"] = parsed[2] |
||||||
|
result["dx"] = parsed[3] |
||||||
|
result["dy"] = parsed[4] |
||||||
|
result["xlonc"] = parsed[5] |
||||||
|
result["truelat1"] = parsed[6] |
||||||
|
result["earth_radius"] = parsed[7] |
||||||
|
|
||||||
|
return result |
||||||
|
|
||||||
|
|
||||||
|
def _parse_record4(data): |
||||||
|
result = {} |
||||||
|
result["is_wind_earth_rel"] = struct.unpack(">i", data) |
||||||
|
|
||||||
|
return result |
||||||
|
|
||||||
|
|
||||||
|
def _parse_record5(data, nx, ny): |
||||||
|
result = {} |
||||||
|
|
||||||
|
size = nx * ny |
||||||
|
fmt = ">{}f".format(size) |
||||||
|
parsed = struct.unpack(fmt, data) |
||||||
|
|
||||||
|
arr = np.array(parsed, dtype=np.float32) |
||||||
|
result["slab"] = arr.reshape((nx, ny), order="F") |
||||||
|
|
||||||
|
return result |
||||||
|
|
||||||
|
|
||||||
|
_PARSE_MAP = {0: _parse_record1, |
||||||
|
1: _parse_record2, |
||||||
|
2: _parse_record3, |
||||||
|
3: _parse_record4, |
||||||
|
4: _parse_record5} |
||||||
|
|
||||||
|
|
||||||
|
def parse_record(record_idx, data, iproj=None, nx=None, ny=None): |
||||||
|
|
||||||
|
if record_idx == 0 or record_idx == 1 or record_idx == 3: |
||||||
|
return _PARSE_MAP[record_idx](data) |
||||||
|
elif record_idx == 2: |
||||||
|
return _PARSE_MAP[record_idx](data, iproj) |
||||||
|
elif record_idx == 4: |
||||||
|
return _PARSE_MAP[record_idx](data, nx, ny) |
||||||
|
|
||||||
|
|
||||||
|
def read_wps(wps_file, field=None): |
||||||
|
with open(wps_file, "rb") as f: |
||||||
|
wps_params = {} |
||||||
|
record_list = [] |
||||||
|
raw_data = f.read() |
||||||
|
|
||||||
|
record_size_idx = 0 |
||||||
|
end_of_file_idx = len(raw_data) - 1 |
||||||
|
|
||||||
|
while True: |
||||||
|
iproj = None |
||||||
|
nx = None |
||||||
|
ny = None |
||||||
|
keep_record = True |
||||||
|
for record_idx in range(5): |
||||||
|
# Each record has the size (in SIZE_BYTES bytes) at the |
||||||
|
# start and end of each record. This might be compiler |
||||||
|
# dependent though, so this might need to be modified. |
||||||
|
# Also, the WPS files are stored big endian. |
||||||
|
|
||||||
|
record_size = struct.unpack( |
||||||
|
">i", |
||||||
|
raw_data[record_size_idx: record_size_idx + SIZE_BYTES]) |
||||||
|
record_start = record_size_idx + SIZE_BYTES |
||||||
|
record_end = record_start + record_size[0] |
||||||
|
record_data = raw_data[record_start:record_end] |
||||||
|
|
||||||
|
parsed_record = parse_record(record_idx, record_data, iproj, |
||||||
|
nx, ny) |
||||||
|
|
||||||
|
try: |
||||||
|
field_name = parsed_record["field"].strip() |
||||||
|
except KeyError: |
||||||
|
pass |
||||||
|
else: |
||||||
|
if field is not None: |
||||||
|
if field_name != field: |
||||||
|
keep_record = False |
||||||
|
|
||||||
|
try: |
||||||
|
iproj = parsed_record["iproj"] |
||||||
|
except KeyError: |
||||||
|
pass |
||||||
|
|
||||||
|
try: |
||||||
|
nx = parsed_record["nx"] |
||||||
|
except KeyError: |
||||||
|
pass |
||||||
|
|
||||||
|
try: |
||||||
|
ny = parsed_record["ny"] |
||||||
|
except KeyError: |
||||||
|
pass |
||||||
|
|
||||||
|
wps_params.update(parsed_record) |
||||||
|
|
||||||
|
record_size_idx = record_end + SIZE_BYTES |
||||||
|
|
||||||
|
if keep_record: |
||||||
|
record_list.append(WPSData(**wps_params)) |
||||||
|
|
||||||
|
# Repeat for all record slabs |
||||||
|
if record_end + SIZE_BYTES > end_of_file_idx: |
||||||
|
break |
||||||
|
|
||||||
|
return record_list |
||||||
@ -0,0 +1,92 @@ |
|||||||
|
input_file = addfile("/Users/ladwig/Documents/wrf_files/wrfout_d02_2010-06-13_21:00:00.nc", "r") |
||||||
|
|
||||||
|
z = wrf_user_getvar(input_file, "z", 0) ; grid point height |
||||||
|
p = wrf_user_getvar(input_file, "pressure", 0) ; total pressure |
||||||
|
|
||||||
|
dimsz = dimsizes(z) |
||||||
|
pivot = (/ (dimsz(2)-1)/2, (dimsz(1)-1)/2 /) ; pivot point is center of domain |
||||||
|
|
||||||
|
; For the new cross section routine |
||||||
|
xopt = True |
||||||
|
xopt@use_pivot = True |
||||||
|
xopt@angle = 45.0 |
||||||
|
;xopt@levels = |
||||||
|
;xopt@latlon = |
||||||
|
xopt@file_handle = input_file |
||||||
|
;xopt@timeidx = |
||||||
|
xopt@linecoords = True |
||||||
|
|
||||||
|
ht_vertcross = wrf_user_vertcross(z, p, pivot, xopt) |
||||||
|
|
||||||
|
printVarSummary(ht_vertcross) |
||||||
|
print(min(ht_vertcross@lats)) |
||||||
|
print(min(ht_vertcross@lons)) |
||||||
|
print(max(ht_vertcross@lats)) |
||||||
|
print(max(ht_vertcross@lons)) |
||||||
|
|
||||||
|
|
||||||
|
xopt@use_pivot = False |
||||||
|
xopt@angle = 0.0 |
||||||
|
;xopt@levels = |
||||||
|
xopt@latlon = True |
||||||
|
xopt@file_handle = input_file |
||||||
|
xopt@timeidx = 0 |
||||||
|
xopt@linecoords = True |
||||||
|
|
||||||
|
loc_param = (/-104.3632, 32.8562, -95.15308, 40.06575 /) ; pivot point is center of domain |
||||||
|
ht_vertcross2 = wrf_user_vertcross(z, p, loc_param, xopt) |
||||||
|
|
||||||
|
printVarSummary(ht_vertcross2) |
||||||
|
print(min(ht_vertcross2@lats)) |
||||||
|
print(min(ht_vertcross2@lons)) |
||||||
|
print(max(ht_vertcross2@lats)) |
||||||
|
print(max(ht_vertcross2@lons)) |
||||||
|
|
||||||
|
print(ht_vertcross2@lats(190)) |
||||||
|
print(ht_vertcross2@lons(190)) |
||||||
|
|
||||||
|
xopt@use_pivot = True |
||||||
|
xopt@angle = 45.0 |
||||||
|
;xopt@levels = |
||||||
|
xopt@latlon = True |
||||||
|
xopt@file_handle = input_file |
||||||
|
xopt@timeidx = 0 |
||||||
|
xopt@linecoords = True |
||||||
|
|
||||||
|
loc_param := (/-99.98572, 36.54949 /) ; pivot point is center of domain |
||||||
|
ht_vertcross3 = wrf_user_vertcross(z, p, loc_param, xopt) |
||||||
|
|
||||||
|
printVarSummary(ht_vertcross3) |
||||||
|
print(min(ht_vertcross3@lats)) |
||||||
|
print(min(ht_vertcross3@lons)) |
||||||
|
print(max(ht_vertcross3@lats)) |
||||||
|
print(max(ht_vertcross3@lons)) |
||||||
|
|
||||||
|
|
||||||
|
xopt@use_pivot = True |
||||||
|
xopt@angle = 45.0 |
||||||
|
xopt@levels = (/1000., 850., 700., 500., 250. /) |
||||||
|
xopt@latlon = True |
||||||
|
xopt@file_handle = input_file |
||||||
|
xopt@timeidx = 0 |
||||||
|
xopt@linecoords = True |
||||||
|
|
||||||
|
loc_param := (/-99.98572, 36.54949 /) ; pivot point is center of domain |
||||||
|
ht_vertcross4 = wrf_user_vertcross(z, p, loc_param, xopt) |
||||||
|
|
||||||
|
printVarSummary(ht_vertcross4) |
||||||
|
print(min(ht_vertcross4@lats)) |
||||||
|
print(min(ht_vertcross4@lons)) |
||||||
|
print(max(ht_vertcross4@lats)) |
||||||
|
print(max(ht_vertcross4@lons)) |
||||||
|
|
||||||
|
o = True |
||||||
|
o@returnInt = False |
||||||
|
o@useTime = 0 |
||||||
|
l = wrf_user_ll_to_xy(input_file, -99.98572, 36.54949, o) |
||||||
|
print(l) |
||||||
|
|
||||||
|
|
||||||
|
l1 = wrf_user_xy_to_ll(input_file, l(1), l(0), o) |
||||||
|
print(l1) |
||||||
|
|
||||||
@ -0,0 +1,81 @@ |
|||||||
|
a = addfile("wrfout_d03_2017-04-03_06:00:00_ctrl","r") |
||||||
|
|
||||||
|
time = 0 |
||||||
|
refl_10cm = wrf_user_getvar(a,"REFL_10CM",time) |
||||||
|
z = wrf_user_getvar(a, "z", time) |
||||||
|
lat = wrf_user_getvar(a, "lat", time) |
||||||
|
lon = wrf_user_getvar(a, "lon", time) |
||||||
|
|
||||||
|
; convert the lat/lon to x,y |
||||||
|
start_lat = 20.9 |
||||||
|
start_lon = 92.5 |
||||||
|
end_lat = 29.2 |
||||||
|
end_lon = 92.5 |
||||||
|
|
||||||
|
opt = True |
||||||
|
|
||||||
|
start_ij = wrf_user_ll_to_ij(a, start_lon, start_lat, opt) |
||||||
|
start_ij = start_ij - 1 |
||||||
|
|
||||||
|
end_ij = wrf_user_ll_to_ij(a, end_lon, end_lat, opt) |
||||||
|
end_ij = end_ij - 1 |
||||||
|
|
||||||
|
start_end = (/start_ij(0), start_ij(1), end_ij(0), end_ij(1)/) |
||||||
|
|
||||||
|
lat_line = wrf_user_intrp2d(lat,start_end,0.0,True) |
||||||
|
nlat = dimsizes(lat_line) |
||||||
|
|
||||||
|
lon_line = wrf_user_intrp2d(lon,start_end,0.0,True) |
||||||
|
|
||||||
|
refl_cross = wrf_user_intrp3d(refl_10cm,z,"v",start_end,0.,True) |
||||||
|
|
||||||
|
; Need to make a vertical coordinate by using the same code as the |
||||||
|
; cross section |
||||||
|
|
||||||
|
; Currently, the vertical coordinate is not set, so let's do it |
||||||
|
; manually. This will be fixed in the next version of NCL. |
||||||
|
; If you want to set these levels yourself, you'll need to copy the |
||||||
|
; code I sent before and manually set the levels in the cross section |
||||||
|
; routine, then do it again here. |
||||||
|
|
||||||
|
z_max = max(z) |
||||||
|
z_min = 0. |
||||||
|
dz = 0.01 * z_max |
||||||
|
nlevels = tointeger( z_max/dz ) |
||||||
|
z_var2d = new( (/nlevels/), typeof(z)) |
||||||
|
z_var2d(0) = z_min |
||||||
|
|
||||||
|
do i=1, nlevels-1 |
||||||
|
z_var2d(i) = z_var2d(0)+i*dz |
||||||
|
end do |
||||||
|
|
||||||
|
refl_cross&Vertical = z_var2d |
||||||
|
|
||||||
|
wks = gsn_open_wks("png","cross") |
||||||
|
cmap := read_colormap_file("BlAqGrYeOrReVi200") |
||||||
|
cmap(0,:) = (/0,0,0,0/) ; make first color fully transparent |
||||||
|
|
||||||
|
resx = True |
||||||
|
resx@gsnMaximize = True |
||||||
|
resx@lbLabelAutoStride = True ; default v6.1.0 |
||||||
|
|
||||||
|
resx@cnFillOn = True ; turn on color fill |
||||||
|
resx@cnLinesOn = False ; turn lines on/off ; True is default |
||||||
|
resx@cnLineLabelsOn = False ; turn line labels on/off ; True is default |
||||||
|
resx@cnFillPalette = cmap |
||||||
|
nLabels = 8 ; arbitrary |
||||||
|
resx@tmXBLabels = new(nLabels,"string") |
||||||
|
resx@tmXBMode = "Explicit" |
||||||
|
|
||||||
|
resx@tmXBValues := toint(fspan(0,nlat-1,nLabels)) |
||||||
|
do i=0,nLabels-1 |
||||||
|
x = lon_line(i) |
||||||
|
y = lat_line(i) |
||||||
|
resx@tmXBLabels(i) = sprintf("%5.1f", y)+"~C~"+sprintf("%5.1f", x) |
||||||
|
end do |
||||||
|
|
||||||
|
resx@tiMainString = "Full South-North Grid Line X-Section" |
||||||
|
|
||||||
|
|
||||||
|
plot1 = gsn_csm_contour(wks, refl_cross, resx ) |
||||||
|
|
||||||
@ -0,0 +1,26 @@ |
|||||||
|
|
||||||
|
;---Open file and calculate slp. |
||||||
|
a = addfile("/Users/ladwig/Documents/wrf_files/rotated_pole_test.nc","r") |
||||||
|
|
||||||
|
t2 = wrf_user_getvar(a,"T2",0) |
||||||
|
|
||||||
|
;---Start the graphics |
||||||
|
wks = gsn_open_wks("x11","wrf") |
||||||
|
|
||||||
|
;---Set some graphical resources |
||||||
|
res = True |
||||||
|
res@gsnMaximize = True |
||||||
|
res@cnFillOn = True |
||||||
|
res@tfDoNDCOverlay = True ; This is necessary if you don't |
||||||
|
; set sfXArray/sfYArray |
||||||
|
|
||||||
|
;---Add additional map resources |
||||||
|
res = wrf_map_resources(a,res) |
||||||
|
|
||||||
|
;---Change some of the resources that were set (these were set to "gray") |
||||||
|
res@mpGeophysicalLineColor = "black" |
||||||
|
res@mpNationalLineColor = "black" |
||||||
|
res@mpUSStateLineColor = "black" |
||||||
|
|
||||||
|
plot = gsn_csm_contour_map(wks,t2,res) |
||||||
|
|
||||||
@ -0,0 +1,21 @@ |
|||||||
|
ifils = systemfunc ("ls /Users/ladwig/Documents/wrf_files/wrf_vortex_single/wrfout_d02_*") |
||||||
|
print(ifils) |
||||||
|
a = addfiles(ifils, "r") |
||||||
|
;a = addfile("/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00.nc", "r") |
||||||
|
|
||||||
|
lats := (/22.0, 25.0, 27.0 /) |
||||||
|
lons := (/-90.0, -87.5, -83.75 /) |
||||||
|
|
||||||
|
opt = True |
||||||
|
opt@useTime = -1 |
||||||
|
opt@returnI = False |
||||||
|
xy = wrf_user_ll_to_xy(a, lons, lats, opt) |
||||||
|
|
||||||
|
print(xy) |
||||||
|
|
||||||
|
x_s = (/10, 50, 90 /) |
||||||
|
y_s = (/10, 50, 90 /) |
||||||
|
|
||||||
|
ll = wrf_user_xy_to_ll(a, x_s, y_s, opt) |
||||||
|
print(ll) |
||||||
|
|
||||||
@ -0,0 +1,416 @@ |
|||||||
|
;-------------------------------------------------------------------------------- |
||||||
|
|
||||||
|
undef("wrf_user_ll_to_xy") |
||||||
|
function wrf_user_ll_to_xy( file_handle, longitude:numeric, latitude:numeric, \ |
||||||
|
opts_args:logical ) |
||||||
|
|
||||||
|
; This is the same as wrf_user_ll_to_ij, but returns 0-based indexes |
||||||
|
|
||||||
|
begin |
||||||
|
; |
||||||
|
; As of NCL V6.0.0, wrf_user_ll_to_ij can now handle a file |
||||||
|
; or a list of files. |
||||||
|
; |
||||||
|
if(typeof(file_handle).eq."file") then |
||||||
|
ISFILE = True |
||||||
|
nc_file = file_handle |
||||||
|
elseif(typeof(file_handle).eq."list") then |
||||||
|
ISFILE = False |
||||||
|
nc_file = file_handle[0] |
||||||
|
else |
||||||
|
print("wrf_user_ll_to_xy: Error: the first argument must be a file or a list of files opened with addfile or addfiles") |
||||||
|
return |
||||||
|
end if |
||||||
|
|
||||||
|
opts = opts_args |
||||||
|
useT = get_res_value(opts,"useTime",0) |
||||||
|
returnI= get_res_value(opts,"returnInt",True) |
||||||
|
|
||||||
|
res = True |
||||||
|
res@MAP_PROJ = nc_file@MAP_PROJ |
||||||
|
res@TRUELAT1 = nc_file@TRUELAT1 |
||||||
|
res@TRUELAT2 = nc_file@TRUELAT2 |
||||||
|
res@STAND_LON = nc_file@STAND_LON |
||||||
|
res@DX = nc_file@DX |
||||||
|
res@DY = nc_file@DY |
||||||
|
|
||||||
|
if (res@MAP_PROJ .eq. 6) then |
||||||
|
res@POLE_LAT = nc_file@POLE_LAT |
||||||
|
res@POLE_LON = nc_file@POLE_LON |
||||||
|
res@LATINC = (res@DY*360.)/2.0/3.141592653589793/6370000. |
||||||
|
res@LONINC = (res@DX*360.)/2.0/3.141592653589793/6370000. |
||||||
|
else |
||||||
|
res@POLE_LAT = 90.0 |
||||||
|
res@POLE_LON = 0.0 |
||||||
|
res@LATINC = 0.0 |
||||||
|
res@LONINC = 0.0 |
||||||
|
end if |
||||||
|
|
||||||
|
if(isfilevar(nc_file,"XLAT")) |
||||||
|
if(ISFILE) then |
||||||
|
XLAT = nc_file->XLAT(useT,:,:) |
||||||
|
XLONG = nc_file->XLONG(useT,:,:) |
||||||
|
else |
||||||
|
XLAT = file_handle[:]->XLAT |
||||||
|
XLONG = file_handle[:]->XLONG |
||||||
|
end if |
||||||
|
else |
||||||
|
if(ISFILE) then |
||||||
|
XLAT = nc_file->XLAT_M(useT,:,:) |
||||||
|
XLONG = nc_file->XLONG_M(useT,:,:) |
||||||
|
else |
||||||
|
XLAT = file_handle[:]->XLAT_M |
||||||
|
XLONG = file_handle[:]->XLONG_M |
||||||
|
end if |
||||||
|
end if |
||||||
|
|
||||||
|
|
||||||
|
if(dimsizes(dimsizes(XLAT)).eq.2) then |
||||||
|
; Rank 2 |
||||||
|
res@REF_LAT = XLAT(0,0) |
||||||
|
res@REF_LON = XLONG(0,0) |
||||||
|
else |
||||||
|
; Rank 3 |
||||||
|
res@REF_LAT = XLAT(useT,0,0) |
||||||
|
res@REF_LON = XLONG(useT,0,0) |
||||||
|
end if |
||||||
|
res@KNOWNI = 1.0 |
||||||
|
res@KNOWNJ = 1.0 |
||||||
|
|
||||||
|
loc = wrf_ll_to_ij (longitude, latitude, res) |
||||||
|
loc = loc - 1 |
||||||
|
|
||||||
|
if (dimsizes(dimsizes(loc)) .eq. 1) then |
||||||
|
loc!0 = "x_y" |
||||||
|
elseif (dimsizes(dimsizes(loc)) .eq. 2) then |
||||||
|
loc!0 = "x_y" |
||||||
|
loc!1 = "idx" |
||||||
|
else ; Not currently supported |
||||||
|
loc!0 = "x_y" |
||||||
|
loc!1 = "domain_idx" |
||||||
|
loc!2 = "idx" |
||||||
|
end if |
||||||
|
|
||||||
|
if ( returnI ) then |
||||||
|
loci = new(dimsizes(loc),integer) |
||||||
|
;loci@_FillValue = default_fillvalue("integer") ; was -999 |
||||||
|
loci = tointeger(loc + .5) |
||||||
|
loci!0 = loc!0 |
||||||
|
return(loci) |
||||||
|
else |
||||||
|
return(loc) |
||||||
|
end if |
||||||
|
|
||||||
|
|
||||||
|
end |
||||||
|
|
||||||
|
;-------------------------------------------------------------------------------- |
||||||
|
|
||||||
|
undef("wrf_user_xy_to_ll") |
||||||
|
function wrf_user_xy_to_ll( file_handle, x:numeric, y:numeric, \ |
||||||
|
opts_args:logical ) |
||||||
|
|
||||||
|
begin |
||||||
|
; |
||||||
|
; As of NCL V6.0.0, wrf_user_ll_to_ij can now handle a file |
||||||
|
; or a list of files. |
||||||
|
; |
||||||
|
if(typeof(file_handle).eq."file") then |
||||||
|
ISFILE = True |
||||||
|
nc_file = file_handle |
||||||
|
elseif(typeof(file_handle).eq."list") then |
||||||
|
ISFILE = False |
||||||
|
nc_file = file_handle[0] |
||||||
|
else |
||||||
|
print("wrf_user_xy_to_ll: Error: the first argument must be a file or a list of files opened with addfile or addfiles") |
||||||
|
return |
||||||
|
end if |
||||||
|
|
||||||
|
opts = opts_args |
||||||
|
useT = get_res_value(opts,"useTime",0) |
||||||
|
|
||||||
|
res = True |
||||||
|
res@MAP_PROJ = nc_file@MAP_PROJ |
||||||
|
res@TRUELAT1 = nc_file@TRUELAT1 |
||||||
|
res@TRUELAT2 = nc_file@TRUELAT2 |
||||||
|
res@STAND_LON = nc_file@STAND_LON |
||||||
|
res@DX = nc_file@DX |
||||||
|
res@DY = nc_file@DY |
||||||
|
|
||||||
|
if (res@MAP_PROJ .eq. 6) then |
||||||
|
res@POLE_LAT = nc_file@POLE_LAT |
||||||
|
res@POLE_LON = nc_file@POLE_LON |
||||||
|
res@LATINC = (res@DY*360.)/2.0/3.141592653589793/6370000. |
||||||
|
res@LONINC = (res@DX*360.)/2.0/3.141592653589793/6370000. |
||||||
|
else |
||||||
|
res@POLE_LAT = 90.0 |
||||||
|
res@POLE_LON = 0.0 |
||||||
|
res@LATINC = 0.0 |
||||||
|
res@LONINC = 0.0 |
||||||
|
end if |
||||||
|
|
||||||
|
|
||||||
|
if(isfilevar(nc_file,"XLAT")) then |
||||||
|
if(ISFILE) then |
||||||
|
XLAT = nc_file->XLAT(useT,:,:) |
||||||
|
XLONG = nc_file->XLONG(useT,:,:) |
||||||
|
else |
||||||
|
XLAT = file_handle[:]->XLAT |
||||||
|
XLONG = file_handle[:]->XLONG |
||||||
|
end if |
||||||
|
else |
||||||
|
if(ISFILE) then |
||||||
|
XLAT = nc_file->XLAT_M(useT,:,:) |
||||||
|
XLONG = nc_file->XLONG_M(useT,:,:) |
||||||
|
else |
||||||
|
XLAT = file_handle[:]->XLAT_M |
||||||
|
XLONG = file_handle[:]->XLONG_M |
||||||
|
end if |
||||||
|
end if |
||||||
|
|
||||||
|
if(dimsizes(dimsizes(XLAT)).eq.2) then |
||||||
|
; Rank 2 |
||||||
|
res@REF_LAT = XLAT(0,0) |
||||||
|
res@REF_LON = XLONG(0,0) |
||||||
|
else |
||||||
|
; Rank 3 |
||||||
|
res@REF_LAT = XLAT(useT,0,0) |
||||||
|
res@REF_LON = XLONG(useT,0,0) |
||||||
|
end if |
||||||
|
res@KNOWNI = 1.0 |
||||||
|
res@KNOWNJ = 1.0 |
||||||
|
|
||||||
|
; Convert to 1-based indexes for Fortran |
||||||
|
new_x = x + 1 |
||||||
|
new_y = y + 1 |
||||||
|
|
||||||
|
loc = wrf_ij_to_ll (new_x,new_y,res) |
||||||
|
|
||||||
|
if (dimsizes(dimsizes(loc)) .eq. 1) then |
||||||
|
loc!0 = "lon_lat" |
||||||
|
elseif (dimsizes(dimsizes(loc)) .eq. 2) then |
||||||
|
loc!0 = "lon_lat" |
||||||
|
loc!1 = "idx" |
||||||
|
else ; Not currently supported |
||||||
|
loc!0 = "lon_lat" |
||||||
|
loc!1 = "domain_idx" |
||||||
|
loc!2 = "idx" |
||||||
|
end if |
||||||
|
|
||||||
|
return(loc) |
||||||
|
|
||||||
|
|
||||||
|
end |
||||||
|
|
||||||
|
;-------------------------------------------------------------------------------- |
||||||
|
|
||||||
|
undef("wrf_user_vertcross") |
||||||
|
function wrf_user_vertcross(var3d:numeric, z_in:numeric, \ |
||||||
|
loc_param:numeric, opts:logical ) |
||||||
|
|
||||||
|
; var3d - 3d field to interpolate (all input fields must be unstaggered) |
||||||
|
; z_in - interpolate to this field (either p/z) |
||||||
|
; loc_param - an array of 4 values representing the start point and end point |
||||||
|
; for the cross section (start_x, start_y, end_x, end_y) OR a single |
||||||
|
; point when opt@use_pivot is True representing the pivot point. |
||||||
|
; The values can be in grid coordinates or lat/lon coordinates |
||||||
|
; (start_x = start_lon, start_y = start_lat, ...). If using |
||||||
|
; lat/lon coordinates, then opt@latlon must be True. |
||||||
|
; opts - optional arguments |
||||||
|
; use_pivot - set to True to indicate that loc_param and angle are used, |
||||||
|
; otherwise loc_param is set to 4 values to indicate a start and |
||||||
|
; end point |
||||||
|
; angle - an angle for vertical plots - 90 represent a WE cross section, |
||||||
|
; ignored if use_pivot is False. |
||||||
|
; levels - the vertical levels to use in the same units as z_in. Set to |
||||||
|
; False to automatically generate the number of levels specified |
||||||
|
; by autolevels. |
||||||
|
; latlon - set to True if the values in loc_param are latitude and longitude |
||||||
|
; values rather than grid values |
||||||
|
; file_handle - must be set to a file handle when latlon is True or |
||||||
|
; linecoords is True, otherwise this is ignored. |
||||||
|
; timeidx - the time index to use for moving nests when latlon is True. Set |
||||||
|
; to 0 if the nest is not moving. |
||||||
|
; linecoords - set to True to include the latitude and longitude coordinates |
||||||
|
; for the cross section line in the output attributes. |
||||||
|
; autolevels - set to the desired number of levels when levels are |
||||||
|
; selected automatically (default 100). |
||||||
|
|
||||||
|
begin |
||||||
|
|
||||||
|
use_pivot = get_res_value(opts, "use_pivot", False) |
||||||
|
angle = get_res_value(opts, "angle", 0.0) |
||||||
|
levels = get_res_value(opts, "levels", new(1,integer)) |
||||||
|
latlon = get_res_value(opts, "latlon", False) |
||||||
|
file_handle = get_res_value(opts, "file_handle", 0) |
||||||
|
timeidx = get_res_value(opts, "timeidx", 0) |
||||||
|
linecoords = get_res_value(opts, "linecoords", False) |
||||||
|
nlevels = get_res_value(opts, "autolevels", 100) |
||||||
|
|
||||||
|
dims = dimsizes(var3d) |
||||||
|
nd = dimsizes(dims) |
||||||
|
|
||||||
|
dimX = dims(nd-1) |
||||||
|
dimY = dims(nd-2) |
||||||
|
dimZ = dims(nd-3) |
||||||
|
|
||||||
|
if ( nd .eq. 4 ) then |
||||||
|
z = z_in(0,:,:,:) |
||||||
|
else |
||||||
|
z = z_in |
||||||
|
end if |
||||||
|
|
||||||
|
; Convert latlon to xy coordinates |
||||||
|
|
||||||
|
if (use_pivot) then |
||||||
|
if (latlon) then |
||||||
|
opt = True |
||||||
|
opt@returnInt = True |
||||||
|
opt@useTime = timeidx |
||||||
|
ij := wrf_user_ll_to_xy(file_handle, loc_param(0), loc_param(1), opt) |
||||||
|
start_x = ij(0) |
||||||
|
start_y = ij(1) |
||||||
|
else |
||||||
|
start_x = loc_param(0) |
||||||
|
start_y = loc_param(1) |
||||||
|
end if |
||||||
|
else |
||||||
|
if (latlon) then |
||||||
|
opt = True |
||||||
|
opt@returnInt = True |
||||||
|
opt@useTime = timeidx |
||||||
|
ij := wrf_user_ll_to_xy(file_handle, (/ loc_param(0), loc_param(2) /), (/ loc_param(1), loc_param(3) /), opt) |
||||||
|
start_x = ij(0,0) |
||||||
|
start_y = ij(1,0) |
||||||
|
end_x = ij(0,1) |
||||||
|
end_y = ij(1,1) |
||||||
|
else |
||||||
|
start_x = loc_param(0) |
||||||
|
start_y = loc_param(1) |
||||||
|
end_x = loc_param(2) |
||||||
|
end_y = loc_param(3) |
||||||
|
end if |
||||||
|
end if |
||||||
|
|
||||||
|
; get the lat/lons along the cross section line if requested |
||||||
|
|
||||||
|
; set the cross section line coordinates if requested |
||||||
|
if (linecoords) then |
||||||
|
|
||||||
|
latname = "XLAT" |
||||||
|
lonname = "XLONG" |
||||||
|
if(.not. isfilevar(file_handle,"XLAT")) then |
||||||
|
if(isfilevar(file_handle,"XLAT_M")) then |
||||||
|
latname = "XLAT_M" |
||||||
|
lonname = "XLONG_M" |
||||||
|
end if |
||||||
|
end if |
||||||
|
|
||||||
|
latvar = _get_wrf_var(file_handle, latname, timeidx) |
||||||
|
lonvar = _get_wrf_var(file_handle, lonname, timeidx) |
||||||
|
|
||||||
|
if (use_pivot) then |
||||||
|
loc := (/start_x, start_y/) |
||||||
|
linelats = wrf_user_intrp2d(latvar, loc, angle, False) |
||||||
|
linelons = wrf_user_intrp2d(lonvar, loc, angle, False) |
||||||
|
else |
||||||
|
loc := (/start_x, start_y, end_x, end_y /) |
||||||
|
linelats = wrf_user_intrp2d(latvar, loc, angle, True) |
||||||
|
linelons = wrf_user_intrp2d(lonvar, loc, angle, True) |
||||||
|
end if |
||||||
|
|
||||||
|
end if |
||||||
|
|
||||||
|
; set vertical cross section |
||||||
|
; Note for wrf_user_set_xy, opt is False when pivot and angle used. |
||||||
|
if (use_pivot) then |
||||||
|
xy = wrf_user_set_xy( z, start_x, start_y, \ ; assumes 0-based indexing in v6.5.0 |
||||||
|
0.0, 0.0, angle, False ) |
||||||
|
|
||||||
|
else |
||||||
|
xy = wrf_user_set_xy( z, start_x, start_y, \ ; assumes 0-based indexing in v6.5.0 |
||||||
|
end_x, end_y, \ |
||||||
|
angle, True) |
||||||
|
|
||||||
|
end if |
||||||
|
xp = dimsizes(xy) |
||||||
|
|
||||||
|
|
||||||
|
; first we interp z |
||||||
|
var2dz = wrf_interp_2d_xy( z, xy) |
||||||
|
|
||||||
|
; interp to constant z grid |
||||||
|
if (all(ismissing(levels))) then |
||||||
|
if(var2dz(0,0) .gt. var2dz(1,0) ) then ; monotonically decreasing coordinate |
||||||
|
z_max = floor(max(z)/10)*10 ; bottom value |
||||||
|
z_min = ceil(min(z)/10)*10 ; top value |
||||||
|
dz = (1.0/nlevels) * (z_max - z_min) |
||||||
|
;nlevels = tointeger( (z_max-z_min)/dz) |
||||||
|
z_var2d = new( (/nlevels/), typeof(z)) |
||||||
|
z_var2d(0) = z_max |
||||||
|
dz = -dz |
||||||
|
else |
||||||
|
z_max = max(z) |
||||||
|
z_min = 0. |
||||||
|
dz = (1.0/nlevels) * z_max |
||||||
|
;nlevels = tointeger( z_max/dz ) |
||||||
|
z_var2d = new( (/nlevels/), typeof(z)) |
||||||
|
z_var2d(0) = z_min |
||||||
|
end if |
||||||
|
|
||||||
|
do i=1, nlevels-1 |
||||||
|
z_var2d(i) = z_var2d(0)+i*dz |
||||||
|
end do |
||||||
|
else |
||||||
|
z_var2d = levels |
||||||
|
nlevels = dimsizes(z_var2d) |
||||||
|
end if |
||||||
|
|
||||||
|
; interp the variable |
||||||
|
if ( dimsizes(dims) .eq. 4 ) then |
||||||
|
var2d = new( (/dims(0), nlevels, xp(0)/), typeof(var2dz)) |
||||||
|
do it = 0,dims(0)-1 |
||||||
|
var2dtmp = wrf_interp_2d_xy( var3d(it,:,:,:), xy) |
||||||
|
do i=0,xp(0)-1 |
||||||
|
var2d(it,:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d) |
||||||
|
end do |
||||||
|
end do |
||||||
|
var2d!0 = var3d!0 |
||||||
|
var2d!1 = "vertical" |
||||||
|
var2d!2 = "cross_line_idx" |
||||||
|
else |
||||||
|
var2d = new( (/nlevels, xp(0)/), typeof(var2dz)) |
||||||
|
var2dtmp = wrf_interp_2d_xy( var3d, xy) |
||||||
|
do i=0,xp(0)-1 |
||||||
|
var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i), z_var2d) |
||||||
|
end do |
||||||
|
var2d!0 = "vertical" |
||||||
|
var2d!1 = "cross_line_idx" |
||||||
|
end if |
||||||
|
|
||||||
|
st_x = tointeger(xy(0,0)) ; + 1 (removed 1-based indexing in 6.5.0) |
||||||
|
st_y = tointeger(xy(0,1)) ; + 1 |
||||||
|
ed_x = tointeger(xy(xp(0)-1,0)) ; + 1 |
||||||
|
ed_y = tointeger(xy(xp(0)-1,1)) ; + 1 |
||||||
|
if (.not. use_pivot) then |
||||||
|
var2d@Orientation = "Cross-Section: (" + \ |
||||||
|
st_x + "," + st_y + ") to (" + \ |
||||||
|
ed_x + "," + ed_y + ")" |
||||||
|
else |
||||||
|
var2d@Orientation = "Cross-Section: (" + \ |
||||||
|
st_x + "," + st_y + ") to (" + \ |
||||||
|
ed_x + "," + ed_y + ") ; center=(" + \ |
||||||
|
start_x + "," + start_y + \ |
||||||
|
") ; angle=" + angle |
||||||
|
end if |
||||||
|
|
||||||
|
if (linecoords) then |
||||||
|
var2d@lats = linelats |
||||||
|
var2d@lons = linelons |
||||||
|
end if |
||||||
|
|
||||||
|
var2d&vertical = z_var2d |
||||||
|
|
||||||
|
return(var2d) |
||||||
|
|
||||||
|
end |
||||||
Loading…
Reference in new issue