[yt-svn] commit/yt: ngoldbaum: Merged in MatthewTurk/yt (pull request #2212)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Jun 28 15:09:23 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/1d8f9705433e/
Changeset: 1d8f9705433e
Branch: yt
User: ngoldbaum
Date: 2016-06-28 22:09:12+00:00
Summary: Merged in MatthewTurk/yt (pull request #2212)
Add internal geographic coordinates
Affected #: 5 files
diff -r 3ca8a739e5296422ea7954ce8d2717908be25dd6 -r 1d8f9705433e9ec690848d8698db85e8903335c9 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -72,7 +72,8 @@
CylindricalCoordinateHandler, \
SphericalCoordinateHandler, \
GeographicCoordinateHandler, \
- SpectralCubeCoordinateHandler
+ SpectralCubeCoordinateHandler, \
+ InternalGeographicCoordinateHandler
# We want to support the movie format in the future.
# When such a thing comes to pass, I'll move all the stuff that is contant up
@@ -491,6 +492,8 @@
cls = SphericalCoordinateHandler
elif self.geometry == "geographic":
cls = GeographicCoordinateHandler
+ elif self.geometry == "internal_geographic":
+ cls = InternalGeographicCoordinateHandler
elif self.geometry == "spectral_cube":
cls = SpectralCubeCoordinateHandler
else:
diff -r 3ca8a739e5296422ea7954ce8d2717908be25dd6 -r 1d8f9705433e9ec690848d8698db85e8903335c9 yt/geometry/coordinates/api.py
--- a/yt/geometry/coordinates/api.py
+++ b/yt/geometry/coordinates/api.py
@@ -23,7 +23,8 @@
from .spherical_coordinates import \
SphericalCoordinateHandler
from .geographic_coordinates import \
- GeographicCoordinateHandler
+ GeographicCoordinateHandler, \
+ InternalGeographicCoordinateHandler
from .spec_cube_coordinates import \
SpectralCubeCoordinateHandler
diff -r 3ca8a739e5296422ea7954ce8d2717908be25dd6 -r 1d8f9705433e9ec690848d8698db85e8903335c9 yt/geometry/coordinates/geographic_coordinates.py
--- a/yt/geometry/coordinates/geographic_coordinates.py
+++ b/yt/geometry/coordinates/geographic_coordinates.py
@@ -23,14 +23,17 @@
pixelize_cylinder, pixelize_aitoff
class GeographicCoordinateHandler(CoordinateHandler):
+ radial_axis = "altitude"
name = "geographic"
- def __init__(self, ds, ordering = ('latitude', 'longitude', 'altitude')):
+ def __init__(self, ds, ordering = None):
+ if not ordering:
+ ordering = ('latitude', 'longitude', self.radial_axis)
super(GeographicCoordinateHandler, self).__init__(ds, ordering)
self.image_units = {}
self.image_units[self.axis_id['latitude']] = (None, None)
self.image_units[self.axis_id['longitude']] = (None, None)
- self.image_units[self.axis_id['altitude']] = ('deg', 'deg')
+ self.image_units[self.axis_id[self.radial_axis]] = ('deg', 'deg')
def setup_fields(self, registry):
# return the fields for r, z, theta
@@ -56,11 +59,11 @@
display_field = False,
units = "")
- f1, f2 = _get_coord_fields(self.axis_id['altitude'])
- registry.add_field(("index", "daltitude"), function = f1,
+ f1, f2 = _get_coord_fields(self.axis_id[self.radial_axis])
+ registry.add_field(("index", "d%s" % (self.radial_axis,)), function = f1,
display_field = False,
units = "code_length")
- registry.add_field(("index", "altitude"), function = f2,
+ registry.add_field(("index", self.radial_axis), function = f2,
display_field = False,
units = "code_length")
@@ -77,10 +80,10 @@
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,
+ def _path_radial_axis(field, data):
+ return data["index", "d%s" % self.radial_axis]
+ registry.add_field(("index", "path_element_%s" % self.radial_axis),
+ function = _path_radial_axis,
units = "code_length")
def _path_latitude(field, data):
# We use r here explicitly
@@ -98,21 +101,6 @@
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:
- 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 _latitude_to_theta(field, data):
# latitude runs from -90 to 90
return (data["latitude"] + 90) * np.pi/180.0
@@ -137,12 +125,47 @@
function = _dlongitude_to_dphi,
units = "")
+ self._setup_radial_fields(registry)
+
+ def _setup_radial_fields(self, registry):
+ # This stays here because we don't want to risk the field detector not
+ # properly getting the data_source, etc.
+ def _altitude_to_radius(field, data):
+ surface_height = data.get_field_parameter("surface_height")
+ if surface_height is None:
+ 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 _retrieve_radial_offset(self, data_source = None):
+ # This returns the factor by which the radial field is multiplied and
+ # the scalar its offset by. Typically the "factor" will *only* be
+ # either 1.0 or -1.0. The order will be factor * r + offset.
+ # Altitude is the radius from the central zone minus the radius of the
+ # surface. Depth to radius is negative value of depth plus the
+ # outermost radius.
+ surface_height = None
+ if data_source is not None:
+ surface_height = data_source.get_field_parameter("surface_height")
+ if surface_height is None:
+ if hasattr(self.ds, "surface_height"):
+ surface_height = self.ds.surface_height
+ else:
+ surface_height = self.ds.quan(0.0, "code_length")
+ return surface_height, 1.0
+
def pixelize(self, dimension, data_source, field, bounds, size,
antialias = True, periodic = True):
if self.axis_name[dimension] in ('latitude', 'longitude'):
return self._cyl_pixelize(data_source, field, bounds, size,
antialias, dimension)
- elif self.axis_name[dimension] == 'altitude':
+ elif self.axis_name[dimension] == self.radial_axis:
return self._ortho_pixelize(data_source, field, bounds, size,
antialias, dimension, periodic)
else:
@@ -164,13 +187,8 @@
def _cyl_pixelize(self, data_source, field, bounds, size, antialias,
dimension):
- 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
+ offset, factor = self._retrieve_radial_offset(data_source)
+ r = factor * data_source['py'] + offset
# Because of the axis-ordering, dimensions 0 and 1 both have r as py
# and the angular coordinate as px. But we need to figure out how to
# convert our coordinate back to an actual angle, based on which
@@ -192,32 +210,28 @@
buff = buff.transpose()
return buff
-
def convert_from_cartesian(self, coord):
raise NotImplementedError
def convert_to_cartesian(self, coord):
- if hasattr(self.ds, "surface_height"):
- surface_height = self.ds.surface_height
- else:
- surface_height = self.ds.quan(0.0, "code_length")
+ offset, factor = self._retrieve_radial_offset()
if isinstance(coord, np.ndarray) and len(coord.shape) > 1:
- alt = self.axis_id['altitude']
+ rad = self.axis_id[self.radial_axis]
lon = self.axis_id['longitude']
lat = self.axis_id['latitude']
- r = coord[:,alt] + surface_height
+ r = factor * coord[:,rad] + offset
theta = coord[:,lon] * np.pi/180
phi = coord[:,lat] * np.pi/180
nc = np.zeros_like(coord)
# r, theta, phi
nc[:,lat] = np.cos(phi) * np.sin(theta)*r
nc[:,lon] = np.sin(phi) * np.sin(theta)*r
- nc[:,alt] = np.cos(theta) * r
+ nc[:,rad] = np.cos(theta) * r
else:
a, b, c = coord
theta = b * np.pi/180
phi = a * np.pi/180
- r = surface_height + c
+ r = factor * c + offset
nc = (np.cos(phi) * np.sin(theta)*r,
np.sin(phi) * np.sin(theta)*r,
np.cos(theta) * r)
@@ -248,7 +262,7 @@
'y / \\sin(\mathrm{latitude})'),
self.axis_id['longitude']:
('R', 'z'),
- self.axis_id['altitude']:
+ self.axis_id[self.radial_axis]:
('longitude', 'latitude')}
for i in list(rv.keys()):
rv[self.axis_name[i]] = rv[i]
@@ -272,14 +286,14 @@
center, display_center = super(
GeographicCoordinateHandler, self).sanitize_center(center, axis)
name = self.axis_name[axis]
- if name == 'altitude':
+ if name == self.radial_axis:
display_center = center
elif name == 'latitude':
display_center = (0.0 * display_center[0],
0.0 * display_center[1],
0.0 * display_center[2])
elif name == 'longitude':
- ri = self.axis_id['altitude']
+ ri = self.axis_id[self.radial_axis]
c = (self.ds.domain_right_edge[ri] +
self.ds.domain_left_edge[ri])/2.0
display_center = [0.0 * display_center[0],
@@ -293,17 +307,108 @@
if width is not None:
width = super(GeographicCoordinateHandler, self).sanitize_width(
axis, width, depth)
- elif name == 'altitude':
- width = [self.ds.domain_width[self.y_axis['altitude']],
- self.ds.domain_width[self.x_axis['altitude']]]
+ elif name == self.radial_axis:
+ rax = self.radial_axis
+ width = [self.ds.domain_width[self.y_axis[rax]],
+ self.ds.domain_width[self.x_axis[rax]]]
elif name == 'latitude':
- ri = self.axis_id['altitude']
+ ri = self.axis_id[self.radial_axis]
# Remember, in spherical coordinates when we cut in theta,
# we create a conic section
width = [2.0*self.ds.domain_width[ri],
2.0*self.ds.domain_width[ri]]
elif name == 'longitude':
- ri = self.axis_id['altitude']
+ ri = self.axis_id[self.radial_axis]
width = [self.ds.domain_width[ri],
2.0*self.ds.domain_width[ri]]
return width
+
+class InternalGeographicCoordinateHandler(GeographicCoordinateHandler):
+ radial_axis = "depth"
+ name = "internal_geographic"
+
+ def _setup_radial_fields(self, registry):
+ # Altitude is the radius from the central zone minus the radius of the
+ # surface.
+ def _depth_to_radius(field, data):
+ outer_radius = data.get_field_parameter("outer_radius")
+ if outer_radius is None:
+ if hasattr(data.ds, "outer_radius"):
+ outer_radius = data.ds.outer_radius
+ else:
+ # Otherwise, we assume that the depth goes to full depth,
+ # so we can look at the domain right edge in depth.
+ rax = self.axis_id[self.radial_axis]
+ outer_radius = data.ds.domain_right_edge[rax]
+ return -1.0 * data["depth"] + outer_radius
+ registry.add_field(("index", "r"),
+ function=_depth_to_radius,
+ units = "code_length")
+ registry.alias(("index", "dr"), ("index", "ddepth"))
+
+ def _retrieve_radial_offset(self, data_source = None):
+ # Depth means switching sign and adding to full radius
+ outer_radius = None
+ if data_source is not None:
+ outer_radius = data_source.get_field_parameter("outer_radius")
+ if outer_radius is None:
+ if hasattr(self.ds, "outer_radius"):
+ outer_radius = self.ds.outer_radius
+ else:
+ # Otherwise, we assume that the depth goes to full depth,
+ # so we can look at the domain right edge in depth.
+ rax = self.axis_id[self.radial_axis]
+ outer_radius = self.ds.domain_right_edge[rax]
+ return outer_radius, -1.0
+
+ _x_pairs = (('latitude', 'longitude'),
+ ('longitude', 'latitude'),
+ ('depth', 'latitude'))
+
+ _y_pairs = (('latitude', 'depth'),
+ ('longitude', 'depth'),
+ ('depth', 'longitude'))
+
+ def sanitize_center(self, center, axis):
+ center, display_center = super(
+ GeographicCoordinateHandler, self).sanitize_center(center, axis)
+ name = self.axis_name[axis]
+ if name == self.radial_axis:
+ display_center = center
+ elif name == 'latitude':
+ display_center = (0.0 * display_center[0],
+ 0.0 * display_center[1],
+ 0.0 * display_center[2])
+ elif name == 'longitude':
+ ri = self.axis_id[self.radial_axis]
+ offset, factor = self._retrieve_radial_offset()
+ outermost = factor * self.ds.domain_left_edge[ri] + offset
+ display_center = [0.0 * display_center[0],
+ 0.0 * display_center[1],
+ 0.0 * display_center[2]]
+ display_center[self.axis_id['latitude']] = outermost/2.0
+ return center, display_center
+
+ def sanitize_width(self, axis, width, depth):
+ name = self.axis_name[axis]
+ if width is not None:
+ width = super(GeographicCoordinateHandler, self).sanitize_width(
+ axis, width, depth)
+ elif name == self.radial_axis:
+ rax = self.radial_axis
+ width = [self.ds.domain_width[self.y_axis[rax]],
+ self.ds.domain_width[self.x_axis[rax]]]
+ elif name == 'latitude':
+ ri = self.axis_id[self.radial_axis]
+ # Remember, in spherical coordinates when we cut in theta,
+ # we create a conic section
+ offset, factor = self._retrieve_radial_offset()
+ outermost = factor * self.ds.domain_left_edge[ri] + offset
+ width = [2.0*outermost, 2.0*outermost]
+ elif name == 'longitude':
+ ri = self.axis_id[self.radial_axis]
+ offset, factor = self._retrieve_radial_offset()
+ outermost = factor * self.ds.domain_left_edge[ri] + offset
+ width = [outermost, 2.0*outermost]
+ return width
+
diff -r 3ca8a739e5296422ea7954ce8d2717908be25dd6 -r 1d8f9705433e9ec690848d8698db85e8903335c9 yt/geometry/coordinates/tests/test_geographic_coordinates.py
--- a/yt/geometry/coordinates/tests/test_geographic_coordinates.py
+++ b/yt/geometry/coordinates/tests/test_geographic_coordinates.py
@@ -54,3 +54,48 @@
# We also want to check that our radius is correct
yield assert_equal, dd["index","r"], \
dd["index","altitude"] + ds.surface_height
+
+def test_internal_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 depth of 1000 maximum, which
+ # means our volume will be that of a shell 1000 wide, starting at r of
+ # outer_radius - 1000.
+ ds = fake_amr_ds(geometry="internal_geographic")
+ ds.outer_radius = ds.quan(5000, "code_length")
+ axes = ["latitude", "longitude", "depth"]
+ for i, axis in enumerate(axes):
+ dd = ds.all_data()
+ fi = ("index", axis)
+ fd = ("index", "d%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.outer_radius - ds.domain_right_edge[2]
+ outer_r = ds.outer_radius
+ 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_depth"], \
+ dd["index", "ddepth"]
+ yield assert_equal, dd["index", "path_element_depth"], \
+ 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"], \
+ -1.0*dd["index","depth"] + ds.outer_radius
diff -r 3ca8a739e5296422ea7954ce8d2717908be25dd6 -r 1d8f9705433e9ec690848d8698db85e8903335c9 yt/testing.py
--- a/yt/testing.py
+++ b/yt/testing.py
@@ -225,6 +225,8 @@
'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
+ 'internal_geographic' :
+ ( (-90.0, -180.0, 0.0), (90.0, 180.0, 1000.0) ), # latlondep
}
def fake_amr_ds(fields = ("Density",), geometry = "cartesian"):
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