Plot accelerator on the map

A. Petrenko (Novosibirsk, 2019)

In [1]:
import numpy as np
import pandas as pd
import geoviews as gv
import geoviews.tile_sources as gts
from cartopy import crs

gv.extension("bokeh")

Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.
In [2]:
!pwd

/c/!Work/BINP/c-tau/Injector/damping_ring
In [3]:
#%output backend='bokeh'
#
#%opts Curve Spread [width=700 height=300 show_grid=True \
#            default_tools=['box_zoom','pan','wheel_zoom','reset']]
In [4]:
from io import StringIO

out = !sdds2stream results/beamline.mag -col=ElementName,s,Profile -pipe=out
DATA = StringIO("\n".join(out))
df_mag = pd.read_csv(DATA, names=['ElementName', 's', 'Profile'], delim_whitespace=True)
In [5]:
df_mag.tail(3)
Out[5]:
ElementName s Profile
5742 QFW1 131.9443 0.0
5743 QFW1 131.8957 0.0
5744 QFW1 131.9443 0.0
In [6]:
df_wall = pd.read_csv('../../drawings_and_maps/results/tunnel.csv')
In [7]:
df_wall.tail(4)
Out[7]:
Unnamed: 0 Z X
99 99 25.7405 -2.0714
100 100 NaN NaN
101 101 9.2595 -11.5484
102 102 30.6565 -2.0834
In [8]:
x0 = 9252311.415718
y0 = 7332694.100875
rotation_angle=36.57

def elegantZX2geoxy(Z, X, rotation_angle=rotation_angle, x0=x0, y0=y0, latitude=54.849):
    angle = rotation_angle*np.pi/180 # Rotation angle (rad)

    scale = np.cos(latitude * np.pi / 180.0)

    x = (Z*np.cos(angle) - X*np.sin(angle))/scale
    y = (Z*np.sin(angle) + X*np.cos(angle))/scale

    return x0 + x, y0 + y
In [9]:
x,y = elegantZX2geoxy(df_wall['Z'].values, df_wall['X'].values)
In [10]:
%opts Path [width=800 height=500 show_grid=True]
%opts Path.Tunnel (color='red', alpha=0.7, line_width=3)
In [11]:
Tunnel = gv.Path((x,y), kdims=['x', 'y'], group='Tunnel', crs = crs.GOOGLE_MERCATOR)
In [12]:
gts.Wikipedia*Tunnel
Out[12]:

Reading xyz.sdds file:

In [13]:
out = !sdds2stream results/xyz.sdds \
    -col=ElementName,s,X,Y,Z,theta,phi,psi -pipe=out
DATA = StringIO("\n".join(out))
df_xyz = pd.read_csv(DATA, names=['ElementName','s','X','Y','Z','theta','phi','psi'],
                   delim_whitespace=True)
In [14]:
df_xyz['X'] = df_xyz['X'] - 12.5
df_xyz['Z'] = df_xyz['Z'] - 1.85
In [15]:
from scipy import interpolate

theta = interpolate.interp1d(
    df_xyz.s.values, df_xyz.theta.values,
    fill_value=(0, 0), bounds_error=False
)

X0 = interpolate.interp1d(
    df_xyz.s.values, df_xyz.X.values,
    fill_value=(0, 0), bounds_error=False
)

Z0 = interpolate.interp1d(
    df_xyz.s.values, df_xyz.Z.values,
    fill_value=(0, 0), bounds_error=False
)
In [16]:
s = df_mag.s.values
nx =  np.cos(theta(s))
nz = -np.sin(theta(s))

Element_width = 0.5 # m

df_mag['X'] = X0(s) + Element_width*df_mag['Profile']*nx
df_mag['Z'] = Z0(s) + Element_width*df_mag['Profile']*nz
In [17]:
df_mag['x'], df_mag['y'] = elegantZX2geoxy(df_mag.Z,df_mag.X)
In [18]:
df_mag.tail(3)
Out[18]:
ElementName s Profile X Z x y
5742 QFW1 131.9443 0.0 -12.5 -1.850024 9.252322e+06 7.332675e+06
5743 QFW1 131.8957 0.0 -12.5 -1.898624 9.252322e+06 7.332675e+06
5744 QFW1 131.9443 0.0 -12.5 -1.850024 9.252322e+06 7.332675e+06
In [19]:
from bokeh.models import HoverTool
hover = HoverTool(tooltips="@ElementName")
In [20]:
%opts Path.mag2D [tools=[hover]] (color='blue', alpha=0.7)
In [21]:
XY_traj = gv.Path(df_mag, kdims=['x','y'], vdims=['ElementName'], crs = crs.GOOGLE_MERCATOR, group='mag2D')
In [22]:
gts.Wikipedia*XY_traj*Tunnel
Out[22]:
In [23]:
gts.EsriImagery*XY_traj*Tunnel
Out[23]: