[yt-svn] commit/yt: 5 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Apr 16 12:55:05 PDT 2014
5 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/42aa7e01e6a1/
Changeset: 42aa7e01e6a1
Branch: yt-3.0
User: MatthewTurk
Date: 2014-04-10 19:27:13
Summary: Adding geographic coordinate handling.
This is becoming somewhat unweildy. The next stage will be refactoring how
centers are computed, which *needs* to be put off for some time.
Affected #: 4 files
diff -r 794b2beb2c48272fed16de10daa07975309b950f -r 42aa7e01e6a135c2370a81c803a19724966abc92 yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -52,6 +52,8 @@
CylindricalCoordinateHandler
from yt.geometry.spherical_coordinates import \
SphericalCoordinateHandler
+from yt.geometry.geographic_coordinates import \
+ GeographicCoordinateHandler
# 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
@@ -369,6 +371,8 @@
self.coordinates = PolarCoordinateHandler(self)
elif self.geometry == "spherical":
self.coordinates = SphericalCoordinateHandler(self)
+ elif self.geometry == "geographic":
+ self.coordinates = GeographicCoordinateHandler(self)
else:
raise YTGeometryNotSupported(self.geometry)
diff -r 794b2beb2c48272fed16de10daa07975309b950f -r 42aa7e01e6a135c2370a81c803a19724966abc92 yt/geometry/geographic_coordinates.py
--- /dev/null
+++ b/yt/geometry/geographic_coordinates.py
@@ -0,0 +1,200 @@
+"""
+Geographic fields
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+from .coordinate_handler import \
+ CoordinateHandler, \
+ _unknown_coord, \
+ _get_coord_fields
+import yt.visualization._MPL as _MPL
+from yt.utilities.lib.misc_utilities import \
+ pixelize_cylinder, pixelize_aitoff
+
+class GeographicCoordinateHandler(CoordinateHandler):
+
+ def __init__(self, pf, ordering = 'latlonalt'):
+ if ordering != 'latlonalt': raise NotImplementedError
+ super(GeographicCoordinateHandler, self).__init__(pf)
+
+ def setup_fields(self, registry):
+ # return the fields for r, z, theta
+ registry.add_field(("index", "dx"), function=_unknown_coord)
+ registry.add_field(("index", "dy"), function=_unknown_coord)
+ registry.add_field(("index", "dz"), function=_unknown_coord)
+ registry.add_field(("index", "x"), function=_unknown_coord)
+ registry.add_field(("index", "y"), function=_unknown_coord)
+ registry.add_field(("index", "z"), function=_unknown_coord)
+ f1, f2 = _get_coord_fields(0, "")
+ registry.add_field(("index", "dlatitude"), function = f1,
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "latitude"), function = f2,
+ display_field = False,
+ units = "")
+
+ f1, f2 = _get_coord_fields(1, "")
+ registry.add_field(("index", "dlongitude"), function = f1,
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "longitude"), function = f2,
+ display_field = False,
+ units = "")
+
+ f1, f2 = _get_coord_fields(2)
+ registry.add_field(("index", "daltitude"), function = f1,
+ display_field = False,
+ units = "code_length")
+ registry.add_field(("index", "altitude"), function = f2,
+ display_field = False,
+ units = "code_length")
+
+ def _SphericalVolume(field, data):
+ # r**2 sin theta dr dtheta dphi
+ # We can use the transformed coordinates here.
+ vol = data["index", "r"]**2.0
+ vol *= data["index", "dr"]
+ vol *= np.sin(data["index", "theta"])
+ vol *= data["index", "dtheta"]
+ vol *= data["index", "dphi"]
+ return vol
+ registry.add_field(("index", "cell_volume"),
+ function=_SphericalVolume,
+ units = "code_length**3")
+
+ # 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.pf, "surface_height", 0.0)
+ 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
+ registry.add_field(("index", "theta"),
+ function = _longitude_to_theta,
+ units = "")
+ def _dlongitude_to_dtheta(field, data):
+ return data["dlongitude"] * np.pi/180.0
+ registry.add_field(("index", "dtheta"),
+ function = _dlongitude_to_dtheta,
+ units = "")
+
+ def _latitude_to_phi(field, data):
+ # latitude runs from -180 to 180.
+ return (data["latitude"] + 90) * np.pi/180.0
+ registry.add_field(("index", "phi"),
+ function = _latitude_to_phi,
+ units = "")
+ def _dlatitude_to_dphi(field, data):
+ return data["dlatitude"] * np.pi/180.0
+ registry.add_field(("index", "dphi"),
+ function = _dlatitude_to_dphi,
+ units = "")
+
+ def pixelize(self, dimension, data_source, field, bounds, size,
+ antialias = True, periodic = True):
+ if dimension in (0, 1):
+ return self._cyl_pixelize(data_source, field, bounds, size,
+ antialias, dimension)
+ elif dimension == 2:
+ return self._ortho_pixelize(data_source, field, bounds, size,
+ antialias, dimension, periodic)
+ else:
+ raise NotImplementedError
+
+ def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
+ dim, periodic):
+ # We should be using fcoords
+ buff = pixelize_aitoff(data_source["theta"], data_source["dtheta"]/2.0,
+ data_source["phi"], data_source["dphi"]/2.0,
+ size, data_source[field], None, None).transpose()
+ return buff
+
+ 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)
+ buff = pixelize_cylinder(data_source['r'],
+ data_source['dr'] / 2.0,
+ 2.0*np.pi - data_source['theta'],
+ data_source['dtheta'] / 2.0, # half-widths
+ size, data_source[field], bounds,
+ input_img = buff)
+ 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
+ return buff
+
+
+ def convert_from_cartesian(self, coord):
+ raise NotImplementedError
+
+ def convert_to_cartesian(self, coord):
+ raise NotImplementedError
+
+ def convert_to_cylindrical(self, coord):
+ raise NotImplementedError
+
+ def convert_from_cylindrical(self, coord):
+ raise NotImplementedError
+
+ def convert_to_spherical(self, coord):
+ raise NotImplementedError
+
+ def convert_from_spherical(self, coord):
+ raise NotImplementedError
+
+ # Despite being mutables, we uses these here to be clear about how these
+ # are generated and to ensure that they are not re-generated unnecessarily
+ axis_name = { 0 : 'latitude', 1 : 'longitude', 2 : 'altitude',
+ 'latitude' : 'latitude',
+ 'longitude' : 'longitude',
+ 'altitude' : 'altitude',
+ 'Latitude' : 'latitude',
+ 'Longitude' : 'longitude',
+ 'Altitude' : 'altitude',
+ 'lat' : 'latitude',
+ 'lon' : 'longitude',
+ 'alt' : 'altitude' }
+
+ axis_id = { 'latitude' : 0, 'longitude' : 1, 'altitude' : 2,
+ 0 : 0, 1 : 1, 2 : 2}
+
+ x_axis = { 'latitude' : 1, 'longitude' : 0, 'altitude' : 0,
+ 0 : 1, 1 : 0, 2 : 0}
+
+ y_axis = { 'latitude' : 2, 'longitude' : 2, 'altitude' : 1,
+ 0 : 2, 1 : 2, 2 : 1}
+
+ @property
+ def period(self):
+ return self.pf.domain_width
+
diff -r 794b2beb2c48272fed16de10daa07975309b950f -r 42aa7e01e6a135c2370a81c803a19724966abc92 yt/geometry/spherical_coordinates.py
--- a/yt/geometry/spherical_coordinates.py
+++ b/yt/geometry/spherical_coordinates.py
@@ -75,6 +75,7 @@
def pixelize(self, dimension, data_source, field, bounds, size,
antialias = True, periodic = True):
+ self.period
if dimension == 0:
return self._ortho_pixelize(data_source, field, bounds, size,
antialias, dimension, periodic)
@@ -102,18 +103,18 @@
dimension):
if dimension == 1:
buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
+ data_source['dr'] / 2.0,
data_source['phi'],
data_source['dphi'] / 2.0, # half-widths
size, data_source[field], bounds)
elif dimension == 2:
buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
+ data_source['dr'] / 2.0,
data_source['theta'],
data_source['dtheta'] / 2.0, # half-widths
size, data_source[field], bounds)
buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
+ data_source['dr'] / 2.0,
2.0*np.pi - data_source['theta'],
data_source['dtheta'] / 2.0, # half-widths
size, data_source[field], bounds,
diff -r 794b2beb2c48272fed16de10daa07975309b950f -r 42aa7e01e6a135c2370a81c803a19724966abc92 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -177,6 +177,15 @@
# Our default width here is the full domain
width = [pf.domain_right_edge[0]*2.0, pf.domain_right_edge[0]*2.0]
center = pf.arr([0.0, 0.0, 0.0], "code_length")
+ elif pf.geometry == "geographic":
+ c_r = ((pf.domain_right_edge + pf.domain_left_edge)/2.0)[2]
+ center = pf.arr([0.0, 0.0, c_r], "code_length")
+ if axis == 2:
+ # latitude slice
+ width = pf.arr([180, 360], "code_length")
+ else:
+ width = [2.0*(pf.domain_right_edge[2] + pf.surface_height),
+ 2.0*(pf.domain_right_edge[2] + pf.surface_height)]
else:
raise NotImplementedError
bounds = (center[x_dict[axis]]-width[0] / 2,
https://bitbucket.org/yt_analysis/yt/commits/3ea5da4638fb/
Changeset: 3ea5da4638fb
Branch: yt-3.0
User: MatthewTurk
Date: 2014-04-10 19:35:15
Summary: Was confused about "theta slicing."
Affected #: 1 file
diff -r 42aa7e01e6a135c2370a81c803a19724966abc92 -r 3ea5da4638fb411cd8baec862d9b668307745629 yt/geometry/spherical_coordinates.py
--- a/yt/geometry/spherical_coordinates.py
+++ b/yt/geometry/spherical_coordinates.py
@@ -113,12 +113,6 @@
data_source['theta'],
data_source['dtheta'] / 2.0, # half-widths
size, data_source[field], bounds)
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'] / 2.0,
- 2.0*np.pi - data_source['theta'],
- data_source['dtheta'] / 2.0, # half-widths
- size, data_source[field], bounds,
- input_img = buff)
else:
raise RuntimeError
return buff
https://bitbucket.org/yt_analysis/yt/commits/1fb252514fe1/
Changeset: 1fb252514fe1
Branch: yt-3.0
User: MatthewTurk
Date: 2014-04-10 20:26:22
Summary: Speed up aitoff considerably by allowing for faster checks.
Affected #: 1 file
diff -r 3ea5da4638fb411cd8baec862d9b668307745629 -r 1fb252514fe12327f1a1ee74b7b8fdd4c066bbfd yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -18,6 +18,7 @@
cimport numpy as np
cimport cython
cimport libc.math as math
+from fp_utils cimport fmin, fmax
cdef extern from "stdlib.h":
# NOTE that size_t might not be int
@@ -529,6 +530,12 @@
return img
+cdef void aitoff_thetaphi_to_xy(np.float64_t theta, np.float64_t phi,
+ np.float64_t *x, np.float64_t *y):
+ cdef np.float64_t z = math.sqrt(1 + math.cos(phi) * math.cos(theta / 2.0))
+ x[0] = math.cos(phi) * math.sin(theta / 2.0) / z
+ y[0] = math.sin(phi) / z
+
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -539,15 +546,21 @@
buff_size,
np.ndarray[np.float64_t, ndim=1] field,
extents, input_img = None):
-
+ # http://paulbourke.net/geometry/transformationprojection/
+ # longitude is -pi to pi
+ # latitude is -pi/2 to pi/2
+ # z^2 = 1 + cos(latitude) cos(longitude/2)
+ # x = cos(latitude) sin(longitude/2) / z
+ # y = sin(latitude) / z
cdef np.ndarray[np.float64_t, ndim=2] img
cdef int i, j, nf, fi
cdef np.float64_t x, y, z, zb
cdef np.float64_t dx, dy, inside
cdef np.float64_t theta1, dtheta1, phi1, dphi1
- cdef np.float64_t theta0, phi0
+ cdef np.float64_t theta0, phi0, theta_p, dtheta_p, phi_p, dphi_p
cdef np.float64_t PI = np.pi
cdef np.float64_t s2 = math.sqrt(2.0)
+ cdef np.float64_t xmax, ymax, xmin, ymin
nf = field.shape[0]
if input_img is None:
@@ -555,31 +568,66 @@
img[:] = np.nan
else:
img = input_img
+ # Okay, here's our strategy. We compute the bounds in x and y, which will
+ # be a rectangle, and then for each x, y position we check to see if it's
+ # within our theta. This will cost *more* computations of the
+ # (x,y)->(theta,phi) calculation, but because we no longer have to search
+ # through the theta, phi arrays, it should be faster.
dx = 2.0 / (img.shape[0] - 1)
dy = 2.0 / (img.shape[1] - 1)
- for i in range(img.shape[0]):
- x = (-1.0 + i*dx)*s2*2.0
- for j in range(img.shape[1]):
- y = (-1.0 + j * dy)*s2
- zb = (x*x/8.0 + y*y/2.0 - 1.0)
- if zb > 0: continue
- z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
- z = z**0.5
- # Longitude
- phi0 = (2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0))) + PI)
- # Latitude
- # We shift it into co-latitude
- theta0 = (math.asin(z*y) + PI/2.0)
- # Now we just need to figure out which pixel contributes.
- # We do not have a fast search.
- for fi in range(nf):
- theta1 = theta[fi]
- dtheta1 = dtheta[fi]
- if not (theta1 - dtheta1 <= theta0 <= theta1 + dtheta1):
+ for fi in range(nf):
+ theta_p = theta[fi] - PI
+ dtheta_p = dtheta[fi]
+ phi_p = phi[fi] - PI/2.0
+ dphi_p = dphi[fi]
+ # Four transformations
+ aitoff_thetaphi_to_xy(theta_p - dtheta_p, phi_p - dphi_p, &x, &y)
+ xmin = x
+ xmax = x
+ ymin = y
+ ymax = y
+ aitoff_thetaphi_to_xy(theta_p - dtheta_p, phi_p + dphi_p, &x, &y)
+ xmin = fmin(xmin, x)
+ xmax = fmax(xmax, x)
+ ymin = fmin(ymin, y)
+ ymax = fmax(ymax, y)
+ aitoff_thetaphi_to_xy(theta_p + dtheta_p, phi_p - dphi_p, &x, &y)
+ xmin = fmin(xmin, x)
+ xmax = fmax(xmax, x)
+ ymin = fmin(ymin, y)
+ ymax = fmax(ymax, y)
+ aitoff_thetaphi_to_xy(theta_p + dtheta_p, phi_p + dphi_p, &x, &y)
+ xmin = fmin(xmin, x)
+ xmax = fmax(xmax, x)
+ ymin = fmin(ymin, y)
+ ymax = fmax(ymax, y)
+ # Now we have the (projected rectangular) bounds.
+ xmin = (xmin + 1) # Get this into normalized image coords
+ xmax = (xmax + 1) # Get this into normalized image coords
+ ymin = (ymin + 1) # Get this into normalized image coords
+ ymax = (ymax + 1) # Get this into normalized image coords
+ x0 = <int> (xmin / dx)
+ x1 = <int> (xmax / dx) + 1
+ y0 = <int> (ymin / dy)
+ y1 = <int> (ymax / dy) + 1
+ for i in range(x0, x1):
+ x = (-1.0 + i*dx)*s2*2.0
+ for j in range(y0, y1):
+ y = (-1.0 + j * dy)*s2
+ zb = (x*x/8.0 + y*y/2.0 - 1.0)
+ if zb > 0: continue
+ z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
+ z = z**0.5
+ # Longitude
+ theta0 = 2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0)))
+ # Latitude
+ # We shift it into co-latitude
+ phi0 = math.asin(z*y)
+ # Now we just need to figure out which pixel contributes.
+ # We do not have a fast search.
+ if not (theta_p - dtheta_p <= theta0 <= theta_p + dtheta_p):
continue
- phi1 = phi[fi]
- dphi1 = dphi[fi]
- if not (phi1 - dphi1 <= phi0 <= phi1 + dphi1):
+ if not (phi_p - dphi_p <= phi0 <= phi_p + dphi_p):
continue
img[i, j] = field[fi]
return img
https://bitbucket.org/yt_analysis/yt/commits/d5fff50cf6ae/
Changeset: d5fff50cf6ae
Branch: yt-3.0
User: MatthewTurk
Date: 2014-04-10 21:00:16
Summary: Fix up a few more geographic things.
Affected #: 2 files
diff -r 1fb252514fe12327f1a1ee74b7b8fdd4c066bbfd -r d5fff50cf6ae6d9e75e4587f3323ca39e6e04996 yt/geometry/geographic_coordinates.py
--- a/yt/geometry/geographic_coordinates.py
+++ b/yt/geometry/geographic_coordinates.py
@@ -99,7 +99,7 @@
units = "")
def _latitude_to_phi(field, data):
- # latitude runs from -180 to 180.
+ # latitude runs from -90 to 90
return (data["latitude"] + 90) * np.pi/180.0
registry.add_field(("index", "phi"),
function = _latitude_to_phi,
@@ -123,10 +123,10 @@
def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
dim, periodic):
- # We should be using fcoords
buff = pixelize_aitoff(data_source["theta"], data_source["dtheta"]/2.0,
data_source["phi"], data_source["dphi"]/2.0,
- size, data_source[field], None, None).transpose()
+ size, data_source[field], None,
+ None).transpose()
return buff
def _cyl_pixelize(self, data_source, field, bounds, size, antialias,
@@ -137,12 +137,6 @@
data_source['theta'],
data_source['dtheta'] / 2.0, # half-widths
size, data_source[field], bounds)
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'] / 2.0,
- 2.0*np.pi - data_source['theta'],
- data_source['dtheta'] / 2.0, # half-widths
- size, data_source[field], bounds,
- input_img = buff)
elif dimension == 1:
buff = pixelize_cylinder(data_source['r'],
data_source['dr'] / 2.0,
diff -r 1fb252514fe12327f1a1ee74b7b8fdd4c066bbfd -r d5fff50cf6ae6d9e75e4587f3323ca39e6e04996 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -182,10 +182,11 @@
center = pf.arr([0.0, 0.0, c_r], "code_length")
if axis == 2:
# latitude slice
- width = pf.arr([180, 360], "code_length")
+ width = pf.arr([360, 180], "code_length")
else:
width = [2.0*(pf.domain_right_edge[2] + pf.surface_height),
2.0*(pf.domain_right_edge[2] + pf.surface_height)]
+ center[2] = 0.0
else:
raise NotImplementedError
bounds = (center[x_dict[axis]]-width[0] / 2,
https://bitbucket.org/yt_analysis/yt/commits/85aed0cda094/
Changeset: 85aed0cda094
Branch: yt-3.0
User: MatthewTurk
Date: 2014-04-16 21:54:56
Summary: Merged in MatthewTurk/yt/yt-3.0 (pull request #808)
Geographic coordinates
Affected #: 5 files
diff -r 3d28859d2cec38a349e7335fc04e62978ffcc1a8 -r 85aed0cda09466c9635c6539e36f8bd972ac0eaa yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -52,6 +52,8 @@
CylindricalCoordinateHandler
from yt.geometry.spherical_coordinates import \
SphericalCoordinateHandler
+from yt.geometry.geographic_coordinates import \
+ GeographicCoordinateHandler
# 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
@@ -365,6 +367,8 @@
self.coordinates = PolarCoordinateHandler(self)
elif self.geometry == "spherical":
self.coordinates = SphericalCoordinateHandler(self)
+ elif self.geometry == "geographic":
+ self.coordinates = GeographicCoordinateHandler(self)
else:
raise YTGeometryNotSupported(self.geometry)
diff -r 3d28859d2cec38a349e7335fc04e62978ffcc1a8 -r 85aed0cda09466c9635c6539e36f8bd972ac0eaa yt/geometry/geographic_coordinates.py
--- /dev/null
+++ b/yt/geometry/geographic_coordinates.py
@@ -0,0 +1,194 @@
+"""
+Geographic fields
+
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+from .coordinate_handler import \
+ CoordinateHandler, \
+ _unknown_coord, \
+ _get_coord_fields
+import yt.visualization._MPL as _MPL
+from yt.utilities.lib.misc_utilities import \
+ pixelize_cylinder, pixelize_aitoff
+
+class GeographicCoordinateHandler(CoordinateHandler):
+
+ def __init__(self, pf, ordering = 'latlonalt'):
+ if ordering != 'latlonalt': raise NotImplementedError
+ super(GeographicCoordinateHandler, self).__init__(pf)
+
+ def setup_fields(self, registry):
+ # return the fields for r, z, theta
+ registry.add_field(("index", "dx"), function=_unknown_coord)
+ registry.add_field(("index", "dy"), function=_unknown_coord)
+ registry.add_field(("index", "dz"), function=_unknown_coord)
+ registry.add_field(("index", "x"), function=_unknown_coord)
+ registry.add_field(("index", "y"), function=_unknown_coord)
+ registry.add_field(("index", "z"), function=_unknown_coord)
+ f1, f2 = _get_coord_fields(0, "")
+ registry.add_field(("index", "dlatitude"), function = f1,
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "latitude"), function = f2,
+ display_field = False,
+ units = "")
+
+ f1, f2 = _get_coord_fields(1, "")
+ registry.add_field(("index", "dlongitude"), function = f1,
+ display_field = False,
+ units = "")
+ registry.add_field(("index", "longitude"), function = f2,
+ display_field = False,
+ units = "")
+
+ f1, f2 = _get_coord_fields(2)
+ registry.add_field(("index", "daltitude"), function = f1,
+ display_field = False,
+ units = "code_length")
+ registry.add_field(("index", "altitude"), function = f2,
+ display_field = False,
+ units = "code_length")
+
+ def _SphericalVolume(field, data):
+ # r**2 sin theta dr dtheta dphi
+ # We can use the transformed coordinates here.
+ vol = data["index", "r"]**2.0
+ vol *= data["index", "dr"]
+ vol *= np.sin(data["index", "theta"])
+ vol *= data["index", "dtheta"]
+ vol *= data["index", "dphi"]
+ return vol
+ registry.add_field(("index", "cell_volume"),
+ function=_SphericalVolume,
+ units = "code_length**3")
+
+ # 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.pf, "surface_height", 0.0)
+ 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
+ registry.add_field(("index", "theta"),
+ function = _longitude_to_theta,
+ units = "")
+ def _dlongitude_to_dtheta(field, data):
+ return data["dlongitude"] * np.pi/180.0
+ registry.add_field(("index", "dtheta"),
+ function = _dlongitude_to_dtheta,
+ units = "")
+
+ def _latitude_to_phi(field, data):
+ # latitude runs from -90 to 90
+ return (data["latitude"] + 90) * np.pi/180.0
+ registry.add_field(("index", "phi"),
+ function = _latitude_to_phi,
+ units = "")
+ def _dlatitude_to_dphi(field, data):
+ return data["dlatitude"] * np.pi/180.0
+ registry.add_field(("index", "dphi"),
+ function = _dlatitude_to_dphi,
+ units = "")
+
+ def pixelize(self, dimension, data_source, field, bounds, size,
+ antialias = True, periodic = True):
+ if dimension in (0, 1):
+ return self._cyl_pixelize(data_source, field, bounds, size,
+ antialias, dimension)
+ elif dimension == 2:
+ return self._ortho_pixelize(data_source, field, bounds, size,
+ antialias, dimension, periodic)
+ else:
+ raise NotImplementedError
+
+ def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
+ dim, periodic):
+ buff = pixelize_aitoff(data_source["theta"], data_source["dtheta"]/2.0,
+ data_source["phi"], data_source["dphi"]/2.0,
+ size, data_source[field], None,
+ None).transpose()
+ return buff
+
+ 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
+ return buff
+
+
+ def convert_from_cartesian(self, coord):
+ raise NotImplementedError
+
+ def convert_to_cartesian(self, coord):
+ raise NotImplementedError
+
+ def convert_to_cylindrical(self, coord):
+ raise NotImplementedError
+
+ def convert_from_cylindrical(self, coord):
+ raise NotImplementedError
+
+ def convert_to_spherical(self, coord):
+ raise NotImplementedError
+
+ def convert_from_spherical(self, coord):
+ raise NotImplementedError
+
+ # Despite being mutables, we uses these here to be clear about how these
+ # are generated and to ensure that they are not re-generated unnecessarily
+ axis_name = { 0 : 'latitude', 1 : 'longitude', 2 : 'altitude',
+ 'latitude' : 'latitude',
+ 'longitude' : 'longitude',
+ 'altitude' : 'altitude',
+ 'Latitude' : 'latitude',
+ 'Longitude' : 'longitude',
+ 'Altitude' : 'altitude',
+ 'lat' : 'latitude',
+ 'lon' : 'longitude',
+ 'alt' : 'altitude' }
+
+ axis_id = { 'latitude' : 0, 'longitude' : 1, 'altitude' : 2,
+ 0 : 0, 1 : 1, 2 : 2}
+
+ x_axis = { 'latitude' : 1, 'longitude' : 0, 'altitude' : 0,
+ 0 : 1, 1 : 0, 2 : 0}
+
+ y_axis = { 'latitude' : 2, 'longitude' : 2, 'altitude' : 1,
+ 0 : 2, 1 : 2, 2 : 1}
+
+ @property
+ def period(self):
+ return self.pf.domain_width
+
diff -r 3d28859d2cec38a349e7335fc04e62978ffcc1a8 -r 85aed0cda09466c9635c6539e36f8bd972ac0eaa yt/geometry/spherical_coordinates.py
--- a/yt/geometry/spherical_coordinates.py
+++ b/yt/geometry/spherical_coordinates.py
@@ -75,6 +75,7 @@
def pixelize(self, dimension, data_source, field, bounds, size,
antialias = True, periodic = True):
+ self.period
if dimension == 0:
return self._ortho_pixelize(data_source, field, bounds, size,
antialias, dimension, periodic)
@@ -102,22 +103,16 @@
dimension):
if dimension == 1:
buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
+ data_source['dr'] / 2.0,
data_source['phi'],
data_source['dphi'] / 2.0, # half-widths
size, data_source[field], bounds)
elif dimension == 2:
buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
+ data_source['dr'] / 2.0,
data_source['theta'],
data_source['dtheta'] / 2.0, # half-widths
size, data_source[field], bounds)
- buff = pixelize_cylinder(data_source['r'],
- data_source['dr'],
- 2.0*np.pi - data_source['theta'],
- data_source['dtheta'] / 2.0, # half-widths
- size, data_source[field], bounds,
- input_img = buff)
else:
raise RuntimeError
return buff
diff -r 3d28859d2cec38a349e7335fc04e62978ffcc1a8 -r 85aed0cda09466c9635c6539e36f8bd972ac0eaa yt/utilities/lib/misc_utilities.pyx
--- a/yt/utilities/lib/misc_utilities.pyx
+++ b/yt/utilities/lib/misc_utilities.pyx
@@ -18,6 +18,7 @@
cimport numpy as np
cimport cython
cimport libc.math as math
+from fp_utils cimport fmin, fmax
cdef extern from "stdlib.h":
# NOTE that size_t might not be int
@@ -529,6 +530,12 @@
return img
+cdef void aitoff_thetaphi_to_xy(np.float64_t theta, np.float64_t phi,
+ np.float64_t *x, np.float64_t *y):
+ cdef np.float64_t z = math.sqrt(1 + math.cos(phi) * math.cos(theta / 2.0))
+ x[0] = math.cos(phi) * math.sin(theta / 2.0) / z
+ y[0] = math.sin(phi) / z
+
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
@@ -539,15 +546,21 @@
buff_size,
np.ndarray[np.float64_t, ndim=1] field,
extents, input_img = None):
-
+ # http://paulbourke.net/geometry/transformationprojection/
+ # longitude is -pi to pi
+ # latitude is -pi/2 to pi/2
+ # z^2 = 1 + cos(latitude) cos(longitude/2)
+ # x = cos(latitude) sin(longitude/2) / z
+ # y = sin(latitude) / z
cdef np.ndarray[np.float64_t, ndim=2] img
cdef int i, j, nf, fi
cdef np.float64_t x, y, z, zb
cdef np.float64_t dx, dy, inside
cdef np.float64_t theta1, dtheta1, phi1, dphi1
- cdef np.float64_t theta0, phi0
+ cdef np.float64_t theta0, phi0, theta_p, dtheta_p, phi_p, dphi_p
cdef np.float64_t PI = np.pi
cdef np.float64_t s2 = math.sqrt(2.0)
+ cdef np.float64_t xmax, ymax, xmin, ymin
nf = field.shape[0]
if input_img is None:
@@ -555,31 +568,66 @@
img[:] = np.nan
else:
img = input_img
+ # Okay, here's our strategy. We compute the bounds in x and y, which will
+ # be a rectangle, and then for each x, y position we check to see if it's
+ # within our theta. This will cost *more* computations of the
+ # (x,y)->(theta,phi) calculation, but because we no longer have to search
+ # through the theta, phi arrays, it should be faster.
dx = 2.0 / (img.shape[0] - 1)
dy = 2.0 / (img.shape[1] - 1)
- for i in range(img.shape[0]):
- x = (-1.0 + i*dx)*s2*2.0
- for j in range(img.shape[1]):
- y = (-1.0 + j * dy)*s2
- zb = (x*x/8.0 + y*y/2.0 - 1.0)
- if zb > 0: continue
- z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
- z = z**0.5
- # Longitude
- phi0 = (2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0))) + PI)
- # Latitude
- # We shift it into co-latitude
- theta0 = (math.asin(z*y) + PI/2.0)
- # Now we just need to figure out which pixel contributes.
- # We do not have a fast search.
- for fi in range(nf):
- theta1 = theta[fi]
- dtheta1 = dtheta[fi]
- if not (theta1 - dtheta1 <= theta0 <= theta1 + dtheta1):
+ for fi in range(nf):
+ theta_p = theta[fi] - PI
+ dtheta_p = dtheta[fi]
+ phi_p = phi[fi] - PI/2.0
+ dphi_p = dphi[fi]
+ # Four transformations
+ aitoff_thetaphi_to_xy(theta_p - dtheta_p, phi_p - dphi_p, &x, &y)
+ xmin = x
+ xmax = x
+ ymin = y
+ ymax = y
+ aitoff_thetaphi_to_xy(theta_p - dtheta_p, phi_p + dphi_p, &x, &y)
+ xmin = fmin(xmin, x)
+ xmax = fmax(xmax, x)
+ ymin = fmin(ymin, y)
+ ymax = fmax(ymax, y)
+ aitoff_thetaphi_to_xy(theta_p + dtheta_p, phi_p - dphi_p, &x, &y)
+ xmin = fmin(xmin, x)
+ xmax = fmax(xmax, x)
+ ymin = fmin(ymin, y)
+ ymax = fmax(ymax, y)
+ aitoff_thetaphi_to_xy(theta_p + dtheta_p, phi_p + dphi_p, &x, &y)
+ xmin = fmin(xmin, x)
+ xmax = fmax(xmax, x)
+ ymin = fmin(ymin, y)
+ ymax = fmax(ymax, y)
+ # Now we have the (projected rectangular) bounds.
+ xmin = (xmin + 1) # Get this into normalized image coords
+ xmax = (xmax + 1) # Get this into normalized image coords
+ ymin = (ymin + 1) # Get this into normalized image coords
+ ymax = (ymax + 1) # Get this into normalized image coords
+ x0 = <int> (xmin / dx)
+ x1 = <int> (xmax / dx) + 1
+ y0 = <int> (ymin / dy)
+ y1 = <int> (ymax / dy) + 1
+ for i in range(x0, x1):
+ x = (-1.0 + i*dx)*s2*2.0
+ for j in range(y0, y1):
+ y = (-1.0 + j * dy)*s2
+ zb = (x*x/8.0 + y*y/2.0 - 1.0)
+ if zb > 0: continue
+ z = (1.0 - (x/4.0)**2.0 - (y/2.0)**2.0)
+ z = z**0.5
+ # Longitude
+ theta0 = 2.0*math.atan(z*x/(2.0 * (2.0*z*z-1.0)))
+ # Latitude
+ # We shift it into co-latitude
+ phi0 = math.asin(z*y)
+ # Now we just need to figure out which pixel contributes.
+ # We do not have a fast search.
+ if not (theta_p - dtheta_p <= theta0 <= theta_p + dtheta_p):
continue
- phi1 = phi[fi]
- dphi1 = dphi[fi]
- if not (phi1 - dphi1 <= phi0 <= phi1 + dphi1):
+ if not (phi_p - dphi_p <= phi0 <= phi_p + dphi_p):
continue
img[i, j] = field[fi]
return img
diff -r 3d28859d2cec38a349e7335fc04e62978ffcc1a8 -r 85aed0cda09466c9635c6539e36f8bd972ac0eaa yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -177,6 +177,16 @@
# Our default width here is the full domain
width = [pf.domain_right_edge[0]*2.0, pf.domain_right_edge[0]*2.0]
center = pf.arr([0.0, 0.0, 0.0], "code_length")
+ elif pf.geometry == "geographic":
+ c_r = ((pf.domain_right_edge + pf.domain_left_edge)/2.0)[2]
+ center = pf.arr([0.0, 0.0, c_r], "code_length")
+ if axis == 2:
+ # latitude slice
+ width = pf.arr([360, 180], "code_length")
+ else:
+ width = [2.0*(pf.domain_right_edge[2] + pf.surface_height),
+ 2.0*(pf.domain_right_edge[2] + pf.surface_height)]
+ center[2] = 0.0
else:
raise NotImplementedError
bounds = (center[x_dict[axis]]-width[0] / 2,
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