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)
In [39]:
lat = np.array([35.2758333, 35.19277778])
lon = np.array([-80.8261111, -80.9358333])
xy = latLon2shrap(lat, lon)
print(xy)
In [33]:
xy
Out[33]:
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]:
In [41]:
shrap2latlon(x,y)
Out[41]:
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]:
In [56]:
xy_coords = latLon2shrap(np.array(coords.lat), np.array(coords.lon))
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]:
In [67]:
xy_coords = latLon2shrap(np.array(coords.lat), np.array(coords.lon))
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]:
In [ ]: