Courtesy to Marc Schleiss for the R code

LatLon2SHRAP <- function(lat,lon){

nlat <- length(lat)
nlon <- length(lon)
if(nlat!=nlon){stop("lat and lon must have same length")}

d2rad <- 0.01745329
stlat <- 60.0
earth <- 6371.221
clon  <- 15.0

lat_rad <- lat*d2rad
lon_rad <- lon*d2rad
stlat_rad <- stlat*d2rad
clon_rad <- clon * d2rad

sfactor <- (1 + sin(stlat_rad))/(1+ sin(lat_rad))
R <- earth*cos(lat_rad)*sfactor
x <- R*cos(lon_rad + clon_rad)
y <- R*sin(lon_rad + clon_rad)

XY <- matrix(NA_real_,nrow=nlat,ncol=2)
XY[,1] <- (4*x/4.7625) + 1604
XY[,2] <- (4*y/4.7625) + 6404
return(XY)

}


In [1]:
import numpy as np
import pandas as pd

In [38]:
def latLon2shrap(lat, lon):
    nlat = len(lat)
    nlon = len(lon)
    if nlat != nlon:
        print('not the same length')
    d2rad = 0.01745329
    stlat = 60.0
    earth = 6371.221
    clon = 15.0
    
    lat_rad = np.deg2rad(lat)
    lon_rad = np.deg2rad(lon)
    stlat_rad = np.deg2rad(stlat)
    clon_rad = np.deg2rad(clon)
    
    sfactor = (1 + np.sin(stlat_rad))/(1+ np.sin(lat_rad))
    R = earth*np.cos(lat_rad)*sfactor
    x = R*np.cos(lon_rad + clon_rad)
    y = R*np.sin(lon_rad + clon_rad)
    
    XY = np.zeros((nlat,2))
    print(nlat)
    print(XY.shape)
    print(x)
    print(XY[:,1])
    XY[:,0] = (4*x/4.7625) + 1604
    XY[:,1] = (4*y/4.7625) + 6404
    return(XY)

In [32]:
lat = np.array([35.2758333, 35.19277778])
lon = np.array([-80.8261111, -80.9358333])
xy = latLon2shrap(lat, lon)


2
(2, 2)
[ 2519.54465234  2513.24727786]
[ 0.  0.]

In [39]:
lat = np.array([35.2758333, 35.19277778])
lon = np.array([-80.8261111, -80.9358333])
xy = latLon2shrap(lat, lon)
print(xy)


2
(2, 2)
[ 2519.54354898  2513.24617185]
[ 0.  0.]
[[ 3720.15206214  1689.60317925]
 [ 3714.86292649  1677.17800897]]

In [33]:
xy


Out[33]:
array([[ 3720.15298884,  1689.60320794],
       [ 3714.86385542,  1677.17803824]])

In [36]:
x = xy[:,0]
y = xy[:,1]

SHRAP2LatLon <- function(x,y){

## Input:
## x,y = vector with S-HRAP coordinates

## Ouput:
## LATLON = matrix with lat/lon coordinates

d2rad <- 0.01745329
rad2d <- 57.29577950
stlat <- 60
rad   <- 6371.221
clon  <- 15

Nx <- length(x)
Ny <- length(y)

if(Nx!=Ny){stop("shrap_x and shrap_y must have same length")}

x <- 4.7625*(x - 1604)/4
y <- 4.7625*(y - 6404)/4
R <- sqrt(x*x + y*y)

LATLON <- matrix(NA_real_,nrow=Nx,ncol=2)
LATLON[,1] <- 90 - 2*(rad2d*atan(R/(rad*(1+sin(stlat*d2rad)))))
LATLON[,2] <- atan2(y,x)*rad2d - clon

return(LATLON)

}


In [ ]:


In [40]:
def shrap2latlon(x,y):
    d2rad = 0.01745329
    rad2d = 57.29577950
    stlat = 60
    rad   = 6371.221
    clon  = 15
    
    Nx = len(x)
    Ny = len(y)
    
    x = 4.7625*(x - 1604)/4
    y = 4.7625*(y - 6404)/4
    R = np.sqrt(x*x + y*y)
    
    LATLON = np.zeros((Nx,2))
    LATLON[:,0] = 90 - 2*(np.rad2deg(np.arctan(R/(rad*(1+np.sin(np.deg2rad(stlat)))))))
    LATLON[:,1] = np.rad2deg(np.arctan2(y,x)) - clon
    
    return LATLON

In [37]:
shrap2latlon(x,y)


Out[37]:
array([[ 35.27582822, -80.82610158],
       [ 35.19277271, -80.93582376]])

In [41]:
shrap2latlon(x,y)


Out[41]:
array([[ 35.2758301 , -80.8261016 ],
       [ 35.1927746 , -80.93582378]])

In [ ]:


In [ ]:


In [ ]:


In [53]:
coords = pd.read_csv('./radar_sent/radar_lat_lon.csv', header=None)
coords.index = coords.loc[:,0]
coords.drop(coords.columns[0], 1, inplace=True)
coords.columns= ['lat', 'lon']

In [54]:
coords.head()


Out[54]:
lat lon
0
0 36.006 -81.297
1 36.002 -81.287
2 35.999 -81.276
3 35.995 -81.266
4 35.991 -81.256

In [56]:
xy_coords = latLon2shrap(np.array(coords.lat), np.array(coords.lon))


19600
(19600, 2)
[ 2434.84496979  2436.02309764  2437.2455429  ...,  2597.91387578
  2599.13044411  2600.34710995]
[ 0.  0.  0. ...,  0.  0.  0.]

In [59]:
new_coords = pd.DataFrame(xy_coords)

In [60]:
new_coords.to_csv('radar_xy.csv', header=None, index=None)

In [64]:
coords = pd.read_csv('./radar_sent/gauge_lat_lon.csv', header=None)

In [66]:
coords.columns= ['lat', 'lon']
coords.head()


Out[66]:
lat lon
0 35.302925 -80.749975
1 35.331667 -80.826389
2 35.275833 -80.826111
3 35.192778 -80.935833
4 35.278333 -80.892778

In [67]:
xy_coords = latLon2shrap(np.array(coords.lat), np.array(coords.lon))


71
(71, 2)
[ 2525.53669249  2516.5096638   2519.54354898  2513.24617185  2512.87641631
  2526.52507043  2537.32360791  2532.56624403  2526.59178736  2532.65905778
  2516.30245137  2524.76303188  2520.30095136  2508.9983299   2521.14185734
  2530.80965144  2534.43049563  2509.66229547  2522.20733252  2527.74109382
  2513.3825341   2513.89631686  2536.07965145  2504.51580267  2508.83225196
  2515.51318671  2523.95061625  2519.56380318  2550.1736429   2540.33964027
  2532.178324    2542.43416992  2528.9405297   2497.86102658  2500.68290068
  2511.34825971  2506.78394515  2512.03298502  2542.48401226  2491.56459134
  2504.13575372  2502.01455016  2510.28029964  2511.26806914  2508.02793104
  2520.28034503  2529.04117177  2531.12620018  2513.50646598  2518.01881609
  2524.67069541  2499.31938678  2514.18880971  2521.43369925  2522.0041828
  2525.26252482  2537.89105402  2521.41158081  2515.73771124  2519.98083834
  2506.34227925  2504.53644005  2519.94416075  2513.06485841  2534.21768117
  2525.67523797  2544.42110717  2543.84677874  2538.43285811  2534.50856502
  2524.75707597]
[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]

In [69]:
new_coords = pd.DataFrame(xy_coords)
new_coords.to_csv('gauge_xy.csv', header=None, index=None)

In [70]:
new_coords.head()


Out[70]:
0 1
0 3725.185673 1695.147902
1 3717.603917 1695.218869
2 3720.152062 1689.603179
3 3714.862926 1677.178009
4 3714.552371 1687.396215

In [ ]: