forked from 3rdparty/wrf-python
Bill Ladwig
6 years ago
42 changed files with 3012 additions and 1753 deletions
@ -0,0 +1,40 @@
@@ -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 @@
@@ -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 @@
@@ -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() |
||||
@ -0,0 +1,29 @@
@@ -0,0 +1,29 @@
|
||||
import matplotlib.pyplot as plt |
||||
from mpl_toolkits.basemap import Basemap |
||||
|
||||
|
||||
def main(): |
||||
bm = Basemap(projection="rotpole", |
||||
o_lat_p=36.0, |
||||
o_lon_p=180.0, |
||||
llcrnrlat=-10.590603, |
||||
urcrnrlat=46.591976, |
||||
llcrnrlon=-139.08585, |
||||
urcrnrlon=22.661009, |
||||
lon_0=-106.0, |
||||
rsphere=6370000, |
||||
resolution='l') |
||||
|
||||
fig = plt.figure(figsize=(8, 8)) |
||||
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8]) |
||||
|
||||
bm.drawcoastlines(linewidth=.5) |
||||
|
||||
print(bm.proj4string) |
||||
|
||||
plt.savefig("basemap_map.png") |
||||
plt.close(fig) |
||||
|
||||
|
||||
if __name__ == "__main__": |
||||
main() |
||||
@ -1,26 +1,25 @@
@@ -1,26 +1,25 @@
|
||||
# Test snippet for f2py |
||||
import numpy as np |
||||
|
||||
import wrf._wrffortran |
||||
errlen = int(wrf._wrffortran.constants.errlen) |
||||
|
||||
|
||||
a = np.ones((3,3,3)) |
||||
b = np.zeros((3,3,3,3)) |
||||
a = np.ones((3, 3, 3)) |
||||
b = np.zeros((3, 3, 3, 3)) |
||||
c = np.zeros(errlen, "c") |
||||
errstat = np.array(0) |
||||
errmsg = np.zeros(errlen, "c") |
||||
|
||||
c[:] = "Test String" |
||||
|
||||
|
||||
for i in xrange(2): |
||||
outview = b[i,:] |
||||
outview = b[i, :] |
||||
outview = outview.T |
||||
q = wrf._wrffortran.testfunc(a,outview,c,errstat=errstat,errmsg=errmsg) |
||||
print errstat |
||||
|
||||
q = wrf._wrffortran.testfunc(a, outview, c, errstat=errstat, |
||||
errmsg=errmsg) |
||||
print(errstat) |
||||
|
||||
|
||||
str_bytes = (bytes(char).decode("utf-8") for char in errmsg[:]) |
||||
print repr(errmsg) |
||||
print "".join(str_bytes).strip() |
||||
print(repr(errmsg)) |
||||
print("".join(str_bytes).strip()) |
||||
@ -0,0 +1,230 @@
@@ -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 @@
@@ -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 @@
@@ -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 @@
@@ -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 @@
@@ -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 @@
@@ -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 |
||||
@ -1,30 +0,0 @@
@@ -1,30 +0,0 @@
|
||||
|
||||
import matplotlib.pyplot as plt |
||||
from mpl_toolkits.basemap import Basemap |
||||
|
||||
def main(): |
||||
bm = Basemap(projection = "rotpole", |
||||
o_lat_p = 36.0, |
||||
o_lon_p = 180.0, |
||||
llcrnrlat = -10.590603, |
||||
urcrnrlat = 46.591976, |
||||
llcrnrlon = -139.08585, |
||||
urcrnrlon = 22.661009, |
||||
lon_0 = -106.0, |
||||
rsphere = 6370000, |
||||
resolution = 'l') |
||||
|
||||
fig = plt.figure(figsize=(8,8)) |
||||
ax = fig.add_axes([0.1,0.1,0.8,0.8]) |
||||
|
||||
bm.drawcoastlines(linewidth=.5) |
||||
|
||||
print bm.proj4string |
||||
|
||||
plt.savefig("basemap_map.png") |
||||
plt.close(fig) |
||||
|
||||
|
||||
|
||||
if __name__ == "__main__": |
||||
main() |
||||
Loading…
Reference in new issue