Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
import logging
import numpy as np
import pandas as pd
import pyproj
import xarray as xr
from scipy.spatial import cKDTree
logger = logging.getLogger(__name__)
url = 'http://thredds.met.no/thredds/dodsC/meps25files/meps_det_extracted_2_5km_latest.nc'
# using the powerful xarray library to connect to the opendap server
ds = xr.open_dataset(url)
lat = ds['latitude'][:]
lon = ds['longitude'][:]
time_values = ds['time'][:].values
dt_time = pd.to_datetime(time_values)
try:
wind_speed = ds['wind_speed'].load()
except Exception, e:
logger.error('Error in reading wind speed data', exc_info=True)
ds.close()
# in reality geo_object_list is list of power tower objects with coordinates
geo_object_list = []
def geo_object_coordinates():
return None
# collecting x,y coordinates for power towers
line_cartesian_positions = [
(x, y) for geo_object in geo_object_list for
(x, y) in geo_object_coordinates(geo_object)
]
target_points = map(lambda x_: [x_[0], x_[1]], line_cartesian_positions)
projection = pyproj.Proj(
"+proj=utm +zone=33 +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
)
# setting up model tree
lon2d, lat2d = np.meshgrid(lon, lat, sparse=True)
utm_grid = projection(np.ravel(lon2d), np.ravel(lat2d))
model_grid = list(list(zip(*utm_grid)))
# finding nearest point in model grid for all power towers
distance, index = cKDTree(model_grid).query(target_points)
y, x = np.unravel_index(index, lon.shape)
wind_speed_array = wind_speed.isel(x=xr.DataArray(x), y=xr.DataArray(y))[:, 0, :].values
# setting up pointers to geo_object_list
start_ids = [0]
start_ids.extend([len(geo_object) for geo_object in geo_object_list])
start_ids = np.cumsum(start_ids)
# creating dictionary of wind speed pandas dataframes
# each dataframe consists of as many column as there is power towers for the line
wind_speed_dict = {}
for i, geo_object in enumerate(geo_object_list):
wind_speed_dict[geo_object.object_id] = pd.DataFrame(
wind_speed_array[:, start_ids[i]:start_ids[i + 1]],
index=dt_time
)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment