[yt-svn] commit/yt: xarthisius: Merged in MatthewTurk/yt (pull request #1307)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Sat Jan 3 21:29:19 PST 2015
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/821d24988cef/
Changeset: 821d24988cef
Branch: yt
User: xarthisius
Date: 2015-01-04 05:29:10+00:00
Summary: Merged in MatthewTurk/yt (pull request #1307)
Explicitly use path length for projections
Affected #: 17 files
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -236,6 +236,7 @@
else:
raise NotImplementedError(self.method)
self._set_center(center)
+ self._projected_units = {}
if data_source is None: data_source = self.ds.all_data()
for k, v in data_source.field_parameters.items():
if k not in self.field_parameters or \
@@ -289,7 +290,6 @@
if communication_system.communicators[-1].size > 1:
for chunk in self.data_source.chunks([], "io", local_only = False):
self._initialize_chunk(chunk, tree)
- # This needs to be parallel_objects-ified
with self.data_source._field_parameter_state(self.field_parameters):
for chunk in parallel_objects(self.data_source.chunks(
chunk_fields, "io", local_only = True)):
@@ -329,7 +329,6 @@
np.divide(nvals, nwvals[:,None], nvals)
# We now convert to half-widths and center-points
data = {}
- #non_nan = ~np.any(np.isnan(nvals), axis=-1)
code_length = self.ds.domain_width.units
data['px'] = self.ds.arr(px, code_length)
data['py'] = self.ds.arr(py, code_length)
@@ -343,30 +342,8 @@
for fi, field in enumerate(fields):
finfo = self.ds._get_field_info(*field)
mylog.debug("Setting field %s", field)
- units = finfo.units
- # add length units to "projected units" if non-weighted
- # integral projection
- if self.weight_field is None and not self._sum_only and \
- self.method == 'integrate':
- # See _handle_chunk where we mandate cm
- if units == '':
- input_units = "cm"
- else:
- input_units = "(%s) * cm" % units
- else:
- input_units = units
- # Don't forget [non_nan] somewhere here.
- self[field] = YTArray(field_data[fi].ravel(),
- input_units=input_units,
- registry=self.ds.unit_registry)
- # convert units if non-weighted integral projection
- if self.weight_field is None and not self._sum_only and \
- self.method == 'integrate':
- u_obj = Unit(units, registry=self.ds.unit_registry)
- if ((u_obj.is_code_unit or self.ds.no_cgs_equiv_length) and
- not u_obj.is_dimensionless) and input_units != units:
- final_unit = "(%s) * code_length" % units
- self[field].convert_to_units(final_unit)
+ input_units = self._projected_units[field].units
+ self[field] = self.ds.arr(field_data[fi].ravel(), input_units)
for i in data.keys(): self[i] = data.pop(i)
mylog.info("Projection completed")
@@ -381,18 +358,34 @@
def _handle_chunk(self, chunk, fields, tree):
if self.method == "mip" or self._sum_only:
- dl = 1.0
+ dl = self.ds.quan(1.0, "")
else:
# This gets explicitly converted to cm
- dl = chunk.fwidth[:, self.axis]
- dl.convert_to_units("cm")
+ ax_name = self.ds.coordinates.axis_name[self.axis]
+ dl = chunk["index", "path_element_%s" % (ax_name)]
+ # This is done for cases where our path element does not have a CGS
+ # equivalent. Once "preferred units" have been implemented, this
+ # will not be necessary at all, as the final conversion will occur
+ # at the display layer.
+ if not dl.units.is_dimensionless:
+ dl.convert_to_units("cm")
v = np.empty((chunk.ires.size, len(fields)), dtype="float64")
- for i in range(len(fields)):
- v[:,i] = chunk[fields[i]] * dl
+ for i, field in enumerate(fields):
+ d = chunk[field] * dl
+ v[:,i] = d
+ self._projected_units[field] = d.uq
if self.weight_field is not None:
w = chunk[self.weight_field]
np.multiply(v, w[:,None], v)
np.multiply(w, dl, w)
+ for field in fields:
+ # Note that this removes the dl integration, which is
+ # *correct*, but we are not fully self-consistently carrying it
+ # all through. So if interrupted, the process will have
+ # incorrect units assigned in the projected units. This should
+ # not be a problem, since the weight division occurs
+ # self-consistently with unitfree arrays.
+ self._projected_units[field] /= dl.uq
else:
w = np.ones(chunk.ires.size, dtype="float64")
icoords = chunk.icoords
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/fields/field_detector.py
--- a/yt/fields/field_detector.py
+++ b/yt/fields/field_detector.py
@@ -190,6 +190,8 @@
self.requested_parameters.append(param)
if param in ['bulk_velocity', 'center', 'normal']:
return self.ds.arr(np.random.random(3) * 1e-2, self.fp_units[param])
+ elif param in ['surface_height']:
+ return self.ds.quan(0.0, 'code_length')
elif param in ['axis']:
return 0
elif param.startswith("cp_"):
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -32,6 +32,9 @@
registry.add_field(("index", "d%s" % ax), function = f1,
display_field = False,
units = "code_length")
+ registry.add_field(("index", "path_element_%s" % ax), function = f1,
+ display_field = False,
+ units = "code_length")
registry.add_field(("index", "%s" % ax), function = f2,
display_field = False,
units = "code_length")
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/cylindrical_coordinates.py
--- a/yt/geometry/coordinates/cylindrical_coordinates.py
+++ b/yt/geometry/coordinates/cylindrical_coordinates.py
@@ -47,7 +47,7 @@
display_field = False,
units = "code_length")
- f1, f2 = _get_coord_fields(1)
+ f1, f2 = _get_coord_fields(self.axis_id['z'])
registry.add_field(("index", "dz"), function = f1,
display_field = False,
units = "code_length")
@@ -55,7 +55,7 @@
display_field = False,
units = "code_length")
- f1, f2 = _get_coord_fields(2, "")
+ f1, f2 = _get_coord_fields(self.axis_id['theta'], "")
registry.add_field(("index", "dtheta"), function = f1,
display_field = False,
units = "")
@@ -72,6 +72,22 @@
function=_CylindricalVolume,
units = "code_length**3")
+ def _path_r(field, data):
+ return data["index", "dr"]
+ registry.add_field(("index", "path_element_r"),
+ function = _path_r,
+ units = "code_length")
+ def _path_theta(field, data):
+ # Note: this already assumes cell-centered
+ return data["index", "r"] * data["index", "dtheta"]
+ registry.add_field(("index", "path_element_theta"),
+ function = _path_theta,
+ units = "code_length")
+ def _path_z(field, data):
+ return data["index", "dz"]
+ registry.add_field(("index", "path_element_z"),
+ function = _path_z,
+ units = "code_length")
def pixelize(self, dimension, data_source, field, bounds, size,
antialias = True, periodic = True):
@@ -101,10 +117,10 @@
return buff
def _cyl_pixelize(self, data_source, field, bounds, size, antialias):
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
- data_source['theta'],
- data_source['dtheta']/2.0, # half-widths
+ buff = pixelize_cylinder(data_source['px'],
+ data_source['pdx'],
+ data_source['py'],
+ data_source['pdy'],
size, data_source[field], bounds)
return buff
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/geographic_coordinates.py
--- a/yt/geometry/coordinates/geographic_coordinates.py
+++ b/yt/geometry/coordinates/geographic_coordinates.py
@@ -74,40 +74,64 @@
function=_SphericalVolume,
units = "code_length**3")
+ def _path_altitude(field, data):
+ return data["index", "daltitude"]
+ registry.add_field(("index", "path_element_altitude"),
+ function = _path_altitude,
+ units = "code_length")
+ def _path_latitude(field, data):
+ # We use r here explicitly
+ return data["index", "r"] * \
+ data["index", "dlatitude"] * np.pi/180.0
+ registry.add_field(("index", "path_element_latitude"),
+ function = _path_latitude,
+ units = "code_length")
+ def _path_longitude(field, data):
+ # We use r here explicitly
+ return data["index", "r"] \
+ * data["index", "dlongitude"] * np.pi/180.0 \
+ * np.sin((data["index", "latitude"] + 90.0) * np.pi/180.0)
+ registry.add_field(("index", "path_element_longitude"),
+ function = _path_longitude,
+ units = "code_length")
+
# Altitude is the radius from the central zone minus the radius of the
# surface.
def _altitude_to_radius(field, data):
surface_height = data.get_field_parameter("surface_height")
if surface_height is None:
- surface_height = getattr(data.ds, "surface_height", 0.0)
+ if hasattr(data.ds, "surface_height"):
+ surface_height = data.ds.surface_height
+ else:
+ surface_height = data.ds.quan(0.0, "code_length")
return data["altitude"] + surface_height
registry.add_field(("index", "r"),
function=_altitude_to_radius,
units = "code_length")
registry.alias(("index", "dr"), ("index", "daltitude"))
- def _longitude_to_theta(field, data):
- # longitude runs from -180 to 180.
- return (data["longitude"] + 180) * np.pi/180.0
+ def _latitude_to_theta(field, data):
+ # latitude runs from -90 to 90
+ return (data["latitude"] + 90) * np.pi/180.0
registry.add_field(("index", "theta"),
- function = _longitude_to_theta,
+ function = _latitude_to_theta,
units = "")
- def _dlongitude_to_dtheta(field, data):
- return data["dlongitude"] * np.pi/180.0
+ def _dlatitude_to_dtheta(field, data):
+ return data["dlatitude"] * np.pi/180.0
registry.add_field(("index", "dtheta"),
- function = _dlongitude_to_dtheta,
+ function = _dlatitude_to_dtheta,
units = "")
- def _latitude_to_phi(field, data):
- # latitude runs from -90 to 90
- return (data["latitude"] + 90) * np.pi/180.0
+ def _longitude_to_phi(field, data):
+ # longitude runs from -180 to 180
+ return (data["longitude"] + 90) * np.pi/180.0
registry.add_field(("index", "phi"),
- function = _latitude_to_phi,
+ function = _longitude_to_phi,
units = "")
- def _dlatitude_to_dphi(field, data):
- return data["dlatitude"] * np.pi/180.0
+ def _dlongitude_to_dphi(field, data):
+ return data["dlongitude"] * np.pi/180.0
registry.add_field(("index", "dphi"),
- function = _dlatitude_to_dphi,
+ function = _dlongitude_to_dphi,
units = "")
def pixelize(self, dimension, data_source, field, bounds, size,
@@ -131,20 +155,19 @@
def _cyl_pixelize(self, data_source, field, bounds, size, antialias,
dimension):
- if dimension == 0:
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'] / 2.0,
- data_source['theta'],
- data_source['dtheta'] / 2.0, # half-widths
- size, data_source[field], bounds)
- elif dimension == 1:
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'] / 2.0,
- data_source['phi'],
- data_source['dphi'] / 2.0, # half-widths
- size, data_source[field], bounds)
- else:
- raise RuntimeError
+ surface_height = data_source.get_field_parameter("surface_height")
+ if surface_height is None:
+ if hasattr(data_source.ds, "surface_height"):
+ surface_height = data_source.ds.surface_height
+ else:
+ surface_height = data_source.ds.quan(0.0, "code_length")
+ r = data_source['py'] + surface_height
+ # Because of the axis-ordering, dimensions 0 and 1 both have r as py
+ # and the angular coordinate as px.
+ buff = pixelize_cylinder(r, data_source['pdy'],
+ data_source['px'],
+ data_source['pdx'], # half-widths
+ size, data_source[field], bounds)
return buff
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/polar_coordinates.py
--- a/yt/geometry/coordinates/polar_coordinates.py
+++ b/yt/geometry/coordinates/polar_coordinates.py
@@ -15,86 +15,22 @@
#-----------------------------------------------------------------------------
import numpy as np
+from yt.units.yt_array import YTArray
from .coordinate_handler import \
CoordinateHandler, \
_unknown_coord, \
_get_coord_fields
+from .cylindrical_coordinates import CylindricalCoordinateHandler
+import yt.visualization._MPL as _MPL
from yt.utilities.lib.misc_utilities import \
pixelize_cylinder
-class PolarCoordinateHandler(CoordinateHandler):
+class PolarCoordinateHandler(CylindricalCoordinateHandler):
def __init__(self, ds, ordering = 'rtz'):
if ordering != 'rtz': raise NotImplementedError
super(PolarCoordinateHandler, self).__init__(ds)
- def setup_fields(self, registry):
- # return the fields for r, z, theta
- registry.add_field("dx", function=_unknown_coord)
- registry.add_field("dy", function=_unknown_coord)
- registry.add_field("x", function=_unknown_coord)
- registry.add_field("y", function=_unknown_coord)
-
- f1, f2 = _get_coord_fields(0)
- registry.add_field(("index", "dr"), function = f1,
- display_field = False,
- units = "code_length")
- registry.add_field(("index", "r"), function = f2,
- display_field = False,
- units = "code_length")
-
- f1, f2 = _get_coord_fields(1, "")
- registry.add_field(("index", "dtheta"), function = f1,
- display_field = False,
- units = "")
- registry.add_field(("index", "theta"), function = f2,
- display_field = False,
- units = "")
-
- f1, f2 = _get_coord_fields(2)
- registry.add_field(("index", "dz"), function = f1,
- display_field = False,
- units = "code_length")
- registry.add_field(("index", "z"), function = f2,
- display_field = False,
- units = "code_length")
-
-
- def _CylindricalVolume(field, data):
- return data["dtheta"] * data["r"] * data["dr"] * data["dz"]
- registry.add_field("CellVolume", function=_CylindricalVolume)
-
- def pixelize(self, dimension, data_source, field, bounds, size, antialias = True):
- ax_name = self.axis_name[dimension]
- if ax_name in ('r', 'theta'):
- return self._ortho_pixelize(data_source, field, bounds, size,
- antialias)
- elif ax_name == "z":
- return self._polar_pixelize(data_source, field, bounds, size,
- antialias)
- else:
- # Pixelizing along a cylindrical surface is a bit tricky
- raise NotImplementedError
-
-
- def _ortho_pixelize(self, data_source, field, bounds, size, antialias):
- buff = _MPL.Pixelize(data_source['px'], data_source['py'],
- data_source['pdx'], data_source['pdy'],
- data_source[field], size[0], size[1],
- bounds, int(antialias),
- True, self.period).transpose()
- return buff
-
- def _polar_pixelize(self, data_source, field, bounds, size, antialias):
- # Out bounds here will *always* be what plot window thinks are x0, x1,
- # y0, y1, but which will actually be rmin, rmax, thetamin, thetamax.
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
- data_source['theta'],
- data_source['dtheta'] / 2.0, # half-widths
- size, data_source[field], bounds)
- return buff
-
axis_name = { 0 : 'r', 1 : 'theta', 2 : 'z',
'r' : 'r', 'theta' : 'theta', 'z' : 'z',
'R' : 'r', 'Theta' : 'theta', 'Z' : 'z'}
@@ -108,23 +44,6 @@
y_axis = { 'r' : 2, 'theta' : 2, 'z' : 1,
0 : 2, 1 : 2, 2 : 1}
- _image_axis_name = None
- @property
- def image_axis_name(self):
- if self._image_axis_name is not None:
- return self._image_axis_name
- # This is the x and y axes labels that get displayed. For
- # non-Cartesian coordinates, we usually want to override these for
- # Cartesian coordinates, since we transform them.
- rv = {0: ('theta', 'z'),
- 1: ('x', 'y'),
- 2: ('r', 'z')}
- for i in rv.keys():
- rv[self.axis_name[i]] = rv[i]
- rv[self.axis_name[i].upper()] = rv[i]
- self._image_axis_name = rv
- return rv
-
def convert_from_cartesian(self, coord):
return cartesian_to_cylindrical(coord)
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/spec_cube_coordinates.py
--- a/yt/geometry/coordinates/spec_cube_coordinates.py
+++ b/yt/geometry/coordinates/spec_cube_coordinates.py
@@ -17,6 +17,8 @@
import numpy as np
from .cartesian_coordinates import \
CartesianCoordinateHandler
+from .coordinate_handler import \
+ _get_coord_fields
class SpectralCubeCoordinateHandler(CartesianCoordinateHandler):
@@ -58,6 +60,41 @@
self.axis_field = {}
self.axis_field[self.ds.spec_axis] = _spec_axis
+ def setup_fields(self, registry):
+ if self.ds.no_cgs_equiv_length == False:
+ return super(self, SpectralCubeCoordinateHandler
+ ).setup_fields(registry)
+ for axi, ax in enumerate(self.axis_name):
+ f1, f2 = _get_coord_fields(axi)
+ def _get_length_func():
+ def _length_func(field, data):
+ # Just use axis 0
+ rv = data.ds.arr(data.fcoords[...,0].copy(), field.units)
+ rv[:] = 1.0
+ return rv
+ return _length_func
+ registry.add_field(("index", "d%s" % ax), function = f1,
+ display_field = False,
+ units = "code_length")
+ registry.add_field(("index", "path_element_%s" % ax),
+ function = _get_length_func(),
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "%s" % ax), function = f2,
+ display_field = False,
+ units = "code_length")
+ def _cell_volume(field, data):
+ rv = data["index", "dx"].copy(order='K')
+ rv *= data["index", "dy"]
+ rv *= data["index", "dz"]
+ return rv
+ registry.add_field(("index", "cell_volume"), function=_cell_volume,
+ display_field=False, units = "code_length**3")
+ registry.check_derived_fields(
+ [("index", "dx"), ("index", "dy"), ("index", "dz"),
+ ("index", "x"), ("index", "y"), ("index", "z"),
+ ("index", "cell_volume")])
+
def convert_to_cylindrical(self, coord):
raise NotImplementedError
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/spherical_coordinates.py
--- a/yt/geometry/coordinates/spherical_coordinates.py
+++ b/yt/geometry/coordinates/spherical_coordinates.py
@@ -73,6 +73,26 @@
function=_SphericalVolume,
units = "code_length**3")
+ def _path_r(field, data):
+ return data["index", "dr"]
+ registry.add_field(("index", "path_element_r"),
+ function = _path_r,
+ units = "code_length")
+ def _path_theta(field, data):
+ # Note: this already assumes cell-centered
+ return data["index", "r"] * data["index", "dtheta"]
+ registry.add_field(("index", "path_element_theta"),
+ function = _path_theta,
+ units = "code_length")
+ def _path_phi(field, data):
+ # Note: this already assumes cell-centered
+ return data["index", "r"] \
+ * data["index", "dphi"] \
+ * np.sin(data["index", "theta"])
+ registry.add_field(("index", "path_element_phi"),
+ function = _path_phi,
+ units = "code_length")
+
def pixelize(self, dimension, data_source, field, bounds, size,
antialias = True, periodic = True):
self.period
@@ -87,31 +107,27 @@
def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
dim, periodic):
- # We should be using fcoords
- period = self.period[:2].copy() # dummy here
- period[0] = self.period[self.x_axis[dim]]
- period[1] = self.period[self.y_axis[dim]]
- period = period.in_units("code_length").d
- buff = _MPL.Pixelize(data_source['px'], data_source['py'],
- data_source['pdx'], data_source['pdy'],
- data_source[field], size[0], size[1],
- bounds, int(antialias),
- period, int(periodic)).transpose()
+ buff = pixelize_aitoff(data_source["py"], data_source["pdy"],
+ data_source["px"], data_source["pdx"],
+ size, data_source[field], None,
+ None, theta_offset = 0,
+ phi_offset = 0).transpose()
return buff
+
def _cyl_pixelize(self, data_source, field, bounds, size, antialias,
dimension):
if dimension == 1:
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'] / 2.0,
- data_source['phi'],
- data_source['dphi'] / 2.0, # half-widths
+ buff = pixelize_cylinder(data_source['px'],
+ data_source['pdx'],
+ data_source['py'],
+ data_source['pdy'],
size, data_source[field], bounds)
elif dimension == 2:
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'] / 2.0,
- data_source['theta'],
- data_source['dtheta'] / 2.0, # half-widths
+ buff = pixelize_cylinder(data_source['px'],
+ data_source['pdx'],
+ data_source['py'],
+ data_source['pdy'],
size, data_source[field], bounds)
buff = buff.transpose()
else:
@@ -219,3 +235,17 @@
width = [self.ds.domain_width[0],
2.0*self.ds.domain_width[0]]
return width
+
+ def _sanity_check(self):
+ """This prints out a handful of diagnostics that help verify the
+ dataset is well formed."""
+ # We just check a few things here.
+ dd = self.ds.all_data()
+ r0 = self.ds.domain_left_edge[0]
+ r1 = self.ds.domain_right_edge[0]
+ v1 = 4.0 * np.pi / 3.0 * (r1**3 - r0**3)
+ print "Total volume should be 4*pi*r**3 = %0.16e" % (v1)
+ v2 = dd.quantities.total_quantity("cell_volume")
+ print "Actual volume is %0.16e" % (v2)
+ print "Relative difference: %0.16e" % (np.abs(v2-v1)/(v2+v1))
+
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/tests/test_cartesian_coordinates.py
--- /dev/null
+++ b/yt/geometry/coordinates/tests/test_cartesian_coordinates.py
@@ -0,0 +1,27 @@
+# Some tests for the Cartesian coordinates handler
+
+from yt.testing import *
+
+# Our canonical tests are that we can access all of our fields and we can
+# compute our volume correctly.
+
+def test_cartesian_coordinates():
+ # We're going to load up a simple AMR grid and check its volume
+ # calculations and path length calculations.
+ ds = fake_amr_ds()
+ axes = list(set(ds.coordinates.axis_name.values()))
+ axes.sort()
+ for i, axis in enumerate(axes):
+ dd = ds.all_data()
+ fi = ("index", axis)
+ fd = ("index", "d%s" % axis)
+ fp = ("index", "path_element_%s" % axis)
+ ma = np.argmax(dd[fi])
+ yield assert_equal, dd[fi][ma] + dd[fd][ma] / 2.0, ds.domain_right_edge[i]
+ mi = np.argmin(dd[fi])
+ yield assert_equal, dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i]
+ yield assert_equal, dd[fd].min(), ds.index.get_smallest_dx()
+ yield assert_equal, dd[fd].max(), (ds.domain_width/ds.domain_dimensions)[i]
+ yield assert_equal, dd[fd], dd[fp]
+ yield assert_equal, dd["cell_volume"].sum(dtype="float64"), \
+ ds.domain_width.prod()
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/tests/test_cylindrical_coordinates.py
--- /dev/null
+++ b/yt/geometry/coordinates/tests/test_cylindrical_coordinates.py
@@ -0,0 +1,28 @@
+# Some tests for the Cylindrical coordinates handler
+
+from yt.testing import *
+
+# Our canonical tests are that we can access all of our fields and we can
+# compute our volume correctly.
+
+def test_cylindrical_coordinates():
+ # We're going to load up a simple AMR grid and check its volume
+ # calculations and path length calculations.
+ ds = fake_amr_ds(geometry="cylindrical")
+ axes = ["r", "z", "theta"]
+ for i, axis in enumerate(axes):
+ dd = ds.all_data()
+ fi = ("index", axis)
+ fd = ("index", "d%s" % axis)
+ fp = ("index", "path_element_%s" % axis)
+ ma = np.argmax(dd[fi])
+ yield assert_equal, dd[fi][ma] + dd[fd][ma] / 2.0, ds.domain_right_edge[i].d
+ mi = np.argmin(dd[fi])
+ yield assert_equal, dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i].d
+ yield assert_equal, dd[fd].max(), (ds.domain_width/ds.domain_dimensions)[i].d
+ yield assert_almost_equal, dd["cell_volume"].sum(dtype="float64"), \
+ np.pi*ds.domain_width[0]**2 * ds.domain_width[1]
+ yield assert_equal, dd["index", "path_element_r"], dd["index", "dr"]
+ yield assert_equal, dd["index", "path_element_z"], dd["index", "dz"]
+ yield assert_equal, dd["index", "path_element_theta"], \
+ dd["index", "r"] * dd["index", "dtheta"]
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/tests/test_geographic_coordinates.py
--- /dev/null
+++ b/yt/geometry/coordinates/tests/test_geographic_coordinates.py
@@ -0,0 +1,52 @@
+# Some tests for the geographic coordinates handler
+
+from yt.testing import *
+
+# Our canonical tests are that we can access all of our fields and we can
+# compute our volume correctly.
+
+def test_geographic_coordinates():
+ # We're going to load up a simple AMR grid and check its volume
+ # calculations and path length calculations.
+
+ # Note that we are setting it up to have an altitude of 1000 maximum, which
+ # means our volume will be that of a shell 1000 wide, starting at r of
+ # whatever our surface_height is set to.
+ ds = fake_amr_ds(geometry="geographic")
+ ds.surface_height = ds.quan(5000, "code_length")
+ axes = ["latitude", "longitude", "altitude"]
+ for i, axis in enumerate(axes):
+ dd = ds.all_data()
+ fi = ("index", axis)
+ fd = ("index", "d%s" % axis)
+ fp = ("index", "path_element_%s" % axis)
+ ma = np.argmax(dd[fi])
+ yield assert_equal, dd[fi][ma] + dd[fd][ma] / 2.0, ds.domain_right_edge[i].d
+ mi = np.argmin(dd[fi])
+ yield assert_equal, dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i].d
+ yield assert_equal, dd[fd].max(), (ds.domain_width/ds.domain_dimensions)[i].d
+ inner_r = ds.surface_height
+ outer_r = ds.surface_height + ds.domain_width[2]
+ yield assert_equal, dd["index","dtheta"], \
+ dd["index","dlatitude"]*np.pi/180.0
+ yield assert_equal, dd["index","dphi"], \
+ dd["index","dlongitude"]*np.pi/180.0
+ # Note our terrible agreement here.
+ yield assert_rel_equal, dd["cell_volume"].sum(dtype="float64"), \
+ (4.0/3.0) * np.pi * (outer_r**3 - inner_r**3), \
+ 3
+ yield assert_equal, dd["index", "path_element_altitude"], \
+ dd["index", "daltitude"]
+ yield assert_equal, dd["index", "path_element_altitude"], \
+ dd["index", "dr"]
+ # Note that latitude corresponds to theta, longitude to phi
+ yield assert_equal, dd["index", "path_element_latitude"], \
+ dd["index", "r"] * \
+ dd["index", "dlatitude"] * np.pi/180.0
+ yield assert_equal, dd["index", "path_element_longitude"], \
+ dd["index", "r"] * \
+ dd["index", "dlongitude"] * np.pi/180.0 * \
+ np.sin((dd["index", "latitude"] + 90.0) * np.pi/180.0)
+ # We also want to check that our radius is correct
+ yield assert_equal, dd["index","r"], \
+ dd["index","altitude"] + ds.surface_height
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/tests/test_polar_coordinates.py
--- /dev/null
+++ b/yt/geometry/coordinates/tests/test_polar_coordinates.py
@@ -0,0 +1,29 @@
+# Some tests for the polar coordinates handler
+# (Pretty similar to cylindrical, but different ordering)
+
+from yt.testing import *
+
+# Our canonical tests are that we can access all of our fields and we can
+# compute our volume correctly.
+
+def test_cylindrical_coordinates():
+ # We're going to load up a simple AMR grid and check its volume
+ # calculations and path length calculations.
+ ds = fake_amr_ds(geometry="polar")
+ axes = ["r", "theta", "z"]
+ for i, axis in enumerate(axes):
+ dd = ds.all_data()
+ fi = ("index", axis)
+ fd = ("index", "d%s" % axis)
+ fp = ("index", "path_element_%s" % axis)
+ ma = np.argmax(dd[fi])
+ yield assert_equal, dd[fi][ma] + dd[fd][ma] / 2.0, ds.domain_right_edge[i].d
+ mi = np.argmin(dd[fi])
+ yield assert_equal, dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i].d
+ yield assert_equal, dd[fd].max(), (ds.domain_width/ds.domain_dimensions)[i].d
+ yield assert_almost_equal, dd["cell_volume"].sum(dtype="float64"), \
+ np.pi*ds.domain_width[0]**2 * ds.domain_width[2]
+ yield assert_equal, dd["index", "path_element_r"], dd["index", "dr"]
+ yield assert_equal, dd["index", "path_element_z"], dd["index", "dz"]
+ yield assert_equal, dd["index", "path_element_theta"], \
+ dd["index", "r"] * dd["index", "dtheta"]
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/geometry/coordinates/tests/test_spherical_coordinates.py
--- /dev/null
+++ b/yt/geometry/coordinates/tests/test_spherical_coordinates.py
@@ -0,0 +1,35 @@
+# Some tests for the Cylindrical coordinates handler
+
+from yt.testing import *
+
+# Our canonical tests are that we can access all of our fields and we can
+# compute our volume correctly.
+
+def test_spherical_coordinates():
+ # We're going to load up a simple AMR grid and check its volume
+ # calculations and path length calculations.
+ ds = fake_amr_ds(geometry="spherical")
+ axes = ["r", "theta", "phi"]
+ for i, axis in enumerate(axes):
+ dd = ds.all_data()
+ fi = ("index", axis)
+ fd = ("index", "d%s" % axis)
+ fp = ("index", "path_element_%s" % axis)
+ ma = np.argmax(dd[fi])
+ yield assert_equal, dd[fi][ma] + dd[fd][ma] / 2.0, ds.domain_right_edge[i].d
+ mi = np.argmin(dd[fi])
+ yield assert_equal, dd[fi][mi] - dd[fd][mi] / 2.0, ds.domain_left_edge[i].d
+ yield assert_equal, dd[fd].max(), (ds.domain_width/ds.domain_dimensions)[i].d
+ # Note that we're using a lot of funny transforms to get to this, so we do
+ # not expect to get actual agreement. This is a bit of a shame, but I
+ # don't think it is avoidable as of right now. Real datasets will almost
+ # certainly be correct, if this is correct to 3 decimel places.
+ yield assert_almost_equal, dd["cell_volume"].sum(dtype="float64"), \
+ (4.0/3.0) * np.pi*ds.domain_width[0]**3, \
+ 3 # Allow for GENEROUS error on the volume ...
+ yield assert_equal, dd["index", "path_element_r"], dd["index", "dr"]
+ yield assert_equal, dd["index", "path_element_theta"], \
+ dd["index", "r"] * dd["index", "dtheta"]
+ yield assert_equal, dd["index", "path_element_phi"], \
+ dd["index", "r"] * dd["index", "dphi"] * \
+ np.sin(dd["index","theta"])
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/testing.py
--- a/yt/testing.py
+++ b/yt/testing.py
@@ -179,11 +179,25 @@
ug = load_uniform_grid(data, ndims, length_unit=length_unit, nprocs=nprocs)
return ug
-def fake_amr_ds(fields = ("Density",)):
+_geom_transforms = {
+ # These are the bounds we want. Cartesian we just assume goes 0 .. 1.
+ 'cartesian' : ( (0.0, 0.0, 0.0), (1.0, 1.0, 1.0) ),
+ 'spherical' : ( (0.0, 0.0, 0.0), (1.0, np.pi, 2*np.pi) ),
+ 'cylindrical': ( (0.0, 0.0, 0.0), (1.0, 1.0, 2.0*np.pi) ), # rzt
+ 'polar' : ( (0.0, 0.0, 0.0), (1.0, 2.0*np.pi, 1.0) ), # rtz
+ 'geographic' : ( (-90.0, -180.0, 0.0), (90.0, 180.0, 1000.0) ), # latlonalt
+}
+
+def fake_amr_ds(fields = ("Density",), geometry = "cartesian"):
from yt.frontends.stream.api import load_amr_grids
+ LE, RE = _geom_transforms[geometry]
+ LE = np.array(LE)
+ RE = np.array(RE)
data = []
for gspec in _amr_grid_index:
level, left_edge, right_edge, dims = gspec
+ left_edge = left_edge * (RE - LE) + LE
+ right_edge = right_edge * (RE - LE) + LE
gdata = dict(level = level,
left_edge = left_edge,
right_edge = right_edge,
@@ -191,7 +205,8 @@
for f in fields:
gdata[f] = np.random.random(dims)
data.append(gdata)
- return load_amr_grids(data, [32, 32, 32], 1.0)
+ bbox = np.array([LE, RE]).T
+ return load_amr_grids(data, [32, 32, 32], 1.0, geometry=geometry, bbox=bbox)
def fake_particle_ds(
fields = ("particle_position_x",
@@ -225,8 +240,6 @@
ds = load_particles(data, 1.0, bbox=bbox)
return ds
-
-
def expand_keywords(keywords, full=False):
"""
expand_keywords is a means for testing all possible keyword
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -579,7 +579,9 @@
np.ndarray[np.float64_t, ndim=1] dphi,
buff_size,
np.ndarray[np.float64_t, ndim=1] field,
- extents, input_img = None):
+ extents, input_img = None,
+ np.float64_t theta_offset = 0.0,
+ np.float64_t phi_offset = 0.0):
# http://paulbourke.net/geometry/transformationprojection/
# longitude is -pi to pi
# latitude is -pi/2 to pi/2
@@ -610,9 +612,9 @@
dx = 2.0 / (img.shape[0] - 1)
dy = 2.0 / (img.shape[1] - 1)
for fi in range(nf):
- theta_p = theta[fi] - PI
+ theta_p = (theta[fi] + theta_offset) - PI
dtheta_p = dtheta[fi]
- phi_p = phi[fi] - PI/2.0
+ phi_p = (phi[fi] + phi_offset) - PI/2.0
dphi_p = dphi[fi]
# Four transformations
aitoff_thetaphi_to_xy(theta_p - dtheta_p, phi_p - dphi_p, &x, &y)
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/visualization/fixed_resolution.py
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -80,7 +80,11 @@
104923.1
"""
_exclude_fields = ('pz','pdz','dx','x','y','z',
- ('index','dx'),('index','x'),('index','y'),('index','z'))
+ 'r', 'dr', 'phi', 'dphi', 'theta', 'dtheta',
+ ('index','dx'),('index','x'),('index','y'),('index','z'),
+ ('index', 'r'), ('index', 'dr'),
+ ('index', 'phi'), ('index', 'dphi'),
+ ('index', 'theta'), ('index', 'dtheta'))
def __init__(self, data_source, bounds, buff_size, antialias = True,
periodic = False):
self.data_source = data_source
diff -r 20d4fcc40283217c698762168f98fbe9eb2e29ec -r 821d24988cef8f9163cd73f8acc11d9910367489 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -80,7 +80,7 @@
from pyparsing import ParseFatalException
def get_window_parameters(axis, center, width, ds):
- if ds.geometry in ("cartesian", "spectral_cube", "spherical"):
+ if ds.geometry in ("cartesian", "spectral_cube"):
width = ds.coordinates.sanitize_width(axis, width, None)
center, display_center = ds.coordinates.sanitize_center(center, axis)
elif ds.geometry in ("polar", "cylindrical"):
@@ -88,6 +88,14 @@
width = [ds.domain_right_edge[0]*2.0, ds.domain_right_edge[0]*2.0]
center = ds.arr([0.0, 0.0, 0.0], "code_length")
display_center = center.copy()
+ elif ds.geometry == "spherical":
+ center, display_center = ds.coordinates.sanitize_center(center, axis)
+ if axis == 0:
+ # latitude slice
+ width = ds.arr([2*np.pi, np.pi], "code_length")
+ else:
+ width = [2.0*ds.domain_right_edge[2], 2.0*ds.domain_right_edge[2]]
+ center[2] = 0.0
elif ds.geometry == "geographic":
c_r = ((ds.domain_right_edge + ds.domain_left_edge)/2.0)[2]
center = ds.arr([0.0, 0.0, c_r], "code_length")
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list