forked from 3rdparty/wrf-python
Bill Ladwig
6 years ago
92 changed files with 14000 additions and 12804 deletions
@ -0,0 +1,63 @@ |
|||||||
|
.. _contrib_guide: |
||||||
|
|
||||||
|
Contributor Guide |
||||||
|
================================= |
||||||
|
|
||||||
|
.. note:: |
||||||
|
|
||||||
|
This contributor guide is written for wrf-python v1.3.x. In the |
||||||
|
not-too-distant future, wrf-python will undergo a significant refactoring |
||||||
|
to remove the wrapt decorators (which don't serialize for dask), but the |
||||||
|
concepts will remain the same as described below. |
||||||
|
|
||||||
|
|
||||||
|
Ways to Contribute |
||||||
|
----------------------------- |
||||||
|
|
||||||
|
Users are encouraged to contribute various ways. This includes: |
||||||
|
|
||||||
|
- Submitting a bug report |
||||||
|
- Submitting bug fixes |
||||||
|
- Submitting new Fortran computational routines |
||||||
|
- Submitting new Python computational routines |
||||||
|
- Submitting fully wrapped computational routines |
||||||
|
|
||||||
|
|
||||||
|
Getting the source code |
||||||
|
------------------------------ |
||||||
|
|
||||||
|
The source code is available on GitHub: |
||||||
|
|
||||||
|
https://github.com/NCAR/wrf-python |
||||||
|
|
||||||
|
To checkout the code:: |
||||||
|
|
||||||
|
git clone https://github.com/NCAR/wrf-python |
||||||
|
|
||||||
|
|
||||||
|
Git Flow |
||||||
|
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
||||||
|
|
||||||
|
This project follows the GitFlow Workflow, which you can read about here if it |
||||||
|
is new to you: |
||||||
|
|
||||||
|
https://leanpub.com/git-flow/read |
||||||
|
|
||||||
|
When you first clone the repository, by default you will be on the 'develop' |
||||||
|
branch, which is what you should use for your development. |
||||||
|
|
||||||
|
|
||||||
|
Pull Requests |
||||||
|
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ |
||||||
|
|
||||||
|
In order to submit changes, you must use GitHub to issue a pull request. |
||||||
|
|
||||||
|
|
||||||
|
Overview of WRF-Python Internals |
||||||
|
---------------------------------- |
||||||
|
|
||||||
|
WRF-Python is a collection of diagnostic and interpolation routines for WRF-ARW |
||||||
|
data. The API consists of a handful of functions |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
@ -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() |
||||||
@ -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 @@ |
|||||||
|
# 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 |
||||||
@ -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