2017-01-04 57 views
0

我正在使用workshop中的一些代碼從最接近我指定座標的座標提取netCDF文件中的數據。當僅使用一組座標,我能提取我需要沒有如下的麻煩值:Python:將列表中的座標傳遞給函數

import numpy as np 
import netCDF4 
from math import pi 
from numpy import cos, sin 

def tunnel_fast(latvar,lonvar,lat0,lon0): 
    ''' 
    Find closest point in a set of (lat,lon) points to specified point 
    latvar - 2D latitude variable from an open netCDF dataset 
    lonvar - 2D longitude variable from an open netCDF dataset 
    lat0,lon0 - query point 
    Returns iy,ix such that the square of the tunnel distance 
    between (latval[it,ix],lonval[iy,ix]) and (lat0,lon0) 
    is minimum. 
    ''' 

    rad_factor = pi/180.0 # for trignometry, need angles in radians 
    # Read latitude and longitude from file into numpy arrays 
    latvals = latvar[:] * rad_factor 
    lonvals = lonvar[:] * rad_factor 
    ny,nx = latvals.shape 
    lat0_rad = lat0 * rad_factor 
    lon0_rad = lon0 * rad_factor 
    # Compute numpy arrays for all values, no loops 
    clat,clon = cos(latvals),cos(lonvals) 
    slat,slon = sin(latvals),sin(lonvals) 
    delX = cos(lat0_rad)*cos(lon0_rad) - clat*clon 
    delY = cos(lat0_rad)*sin(lon0_rad) - clat*slon 
    delZ = sin(lat0_rad) - slat; 
    dist_sq = delX**2 + delY**2 + delZ**2 
    minindex_1d = dist_sq.argmin() # 1D index of minimum element 
    iy_min,ix_min = np.unravel_index(minindex_1d, latvals.shape) 
    return iy_min,ix_min 

ncfile = netCDF4.Dataset('E:\wind_level2_1.nc', 'r') 
latvar = ncfile.variables['latitude'] 
lonvar = ncfile.variables['longitude'] 

#_________GG turbine_________GAD10 Latitude 51.735516, GAD10 Longitude 1.942656 

iy,ix = tunnel_fast(latvar, lonvar, 51.735516, 1.942656) 
print('Closest lat lon:', latvar[iy,ix], lonvar[iy,ix]) 

refLAT=latvar[iy,ix] 
refLON = lonvar[iy,ix] 
#try to find the data for this location 

SARwind = ncfile.variables['sar_wind'][:,:] 
ModelWind = ncfile.variables['model_speed'][:,:] 

print 'iy,ix' #appears to be the index of the value of Lat,lon 
print SARwind[iy,ix] 


ncfile.close() 

現在我通過包含座標coord_list提取套座標的文本文件,試圖循環,查找數據然後移動到列表中的下一組座標。此代碼的工作在它自己的如下:

import csv 
from decimal import Decimal 
with open('Turbine_locs_no_header.csv','rb') as f: 
    reader = csv.reader(f) 
    #coord_list = list(reader) 
    coord_list = [reader] 
    end_row = len(coord_list) 

    lon_ind=1 
    lat_ind=2 

for row in range(0, end_row-1):#end_row - 1 due to the 0 index 
    turbine_lat = coord_list[row][lat_ind] 
    turbine_lon = coord_list[row][lon_ind] 
    turbine_lat = [Decimal(turbine_lat)] 
    print 'lat',turbine_lat, 'lon',turbine_lon, row 

不過,我想從文本文件的原代碼iy,ix = tunnel_fast(latvar, lonvar, 51.94341, 1.922094888)這部分傳遞的是座標,與變量iy, ix = tunnel_fast(latvar, lonvar, turbine_lat, turbine_lon)更換號碼。我試圖通過創建一個功能get_coordinates兩個代碼結合,我收到以下錯誤

File "C:/Users/mm/test_nc_bycoords_GG_turbines_AGW.py", line 65, in <module> 
    get_coordinates(coord_list, latvar, lonvar) 
    File "C:/Users/mm/test_nc_bycoords_GG_turbines_AGW.py", line 51, in get_coordinates 
    iy, ix = tunnel_fast(latvar, lonvar, turbine_lat, turbine_lon) 
    File "C:/Users/mm/test_nc_bycoords_GG_turbines_AGW.py", line 27, in tunnel_fast 
    lat0_rad = lat0 * rad_factor 
TypeError: can't multiply sequence by non-int of type 'float' 

我想這是因爲turbine_latturbine_lon是列表項所以不能使用,但這似乎並沒有被連接到錯誤。無論如何,我知道這段代碼需要更多的工作,但是如果任何人都可以幫助我發現我出錯的地方,那將會非常有幫助。我試圖結合這兩個代碼如下。

import numpy as np 
import netCDF4 
from math import pi 
from numpy import cos, sin 
import csv 

# edited from https://github.com/Unidata/unidata-python-workshop/blob/a56daa50d7b343c7debe93968683613642d6b9f7/notebooks/netcdf-by-coordinates.ipynb 

def tunnel_fast(latvar,lonvar,lat0,lon0): 
    ''' 
    Find closest point in a set of (lat,lon) points to specified point 
    latvar - 2D latitude variable from an open netCDF dataset 
    lonvar - 2D longitude variable from an open netCDF dataset 
    lat0,lon0 - query point 
    Returns iy,ix such that the square of the tunnel distance 
    between (latval[it,ix],lonval[iy,ix]) and (lat0,lon0) 
    is minimum. 
    ''' 

    rad_factor = pi/180.0 # for trignometry, need angles in radians 
    # Read latitude and longitude from file into numpy arrays 
    latvals = latvar[:] * rad_factor 
    lonvals = lonvar[:] * rad_factor 
    ny,nx = latvals.shape 
    lat0_rad = lat0 * rad_factor 
    lon0_rad = lon0 * rad_factor 
    # Compute numpy arrays for all values, no loops 
    clat,clon = cos(latvals),cos(lonvals) 
    slat,slon = sin(latvals),sin(lonvals) 
    delX = cos(lat0_rad)*cos(lon0_rad) - clat*clon 
    delY = cos(lat0_rad)*sin(lon0_rad) - clat*slon 
    delZ = sin(lat0_rad) - slat; 
    dist_sq = delX**2 + delY**2 + delZ**2 
    minindex_1d = dist_sq.argmin() # 1D index of minimum element 
    iy_min,ix_min = np.unravel_index(minindex_1d, latvals.shape) 
    return iy_min,ix_min 
#________________my edits___________________________________________________ 
def get_coordinates(coord_list, latvar, lonvar): 
    "this takes coordinates from a .csv and assigns them to variables" 

    end_row = len(coord_list) 

    lon_ind=1 
    lat_ind=2 

    for row in range(0, end_row-1):#end_row - 1 due to the 0 index 
     turbine_lat = coord_list[row][lat_ind] 
     turbine_lon = coord_list[row][lon_ind] 
     iy, ix = tunnel_fast(latvar, lonvar, turbine_lat, turbine_lon) 
     print('Closest lat lon:', latvar[iy, ix], lonvar[iy, ix]) 

#________________________________________________________________________________________________________________________ 
ncfile = netCDF4.Dataset('NOGAPS_wind_level2_1.nc', 'r') 
latvar = ncfile.variables['latitude'] 
lonvar = ncfile.variables['longitude'] 

#____added in to pass to get coordinates function 
with open('Turbine_locs_no_header.csv','rb') as f: 
    reader = csv.reader(f) 
    coord_list = list(reader) 
#_________take latitude from coordinateas function 

get_coordinates(coord_list, latvar, lonvar) 

#iy,ix = tunnel_fast(latvar, lonvar, turbine_lat, turbine_lon)#get these from the 'assign_coordinates_fromlist.py 
#print('Closest lat lon:', latvar[iy,ix], lonvar[iy,ix]) 

SARwind = ncfile.variables['sar_wind'][:,:] 
ModelWind = ncfile.variables['model_speed'][:,:] 

print 'iy,ix' #appears to be the index of the value of Lat,lon 
print SARwind[iy,ix] 


ncfile.close() 

當我嘗試轉換

+0

您能更新第一個代碼示例中的代碼格式嗎?我會但不是6個角色的更新。 –

+0

@EmettSpeer。我已經在第一部分修復了縮進,還有什麼需要排序? –

+0

我只發現了一個格式問題。 –

回答

0

感謝來自a_guest

這是因爲 <type 'str'>傳遞的lat0lon0一個簡單的問題,以tunnel_fast當它需要<type 'float'>幫助。這似乎來自加載coord_list作爲列表。

with open('Turbine_locs_no_header.csv','rb') as f: 
    reader = csv.reader(f) 
    coord_list = list(reader) 

我所用的解決方法是轉換lat0lon0tunnel_fast

lat0 = float(lat0) 
lon0 = float(lon0) 

開始到彩車我相信有一個更優雅的方式來做到這一點,但它的工作原理。

1

可以解壓使用*args參數列表(見the docs)。在你的情況下,你可以做tunnel_fast(latvar, lonvar, *coord_list[row])。您需要確保coord_list[row]中的參數順序是正確的,如果coord_list[row]包含多於兩個值,則需要適當地對其進行分片。

+0

謝謝了,回溯錯誤與'iy,ix = tunnel_fast(latvar,lonvar,* coord_list [row])'而不是'iy,ix = tunnel_fast(latvar,lonvar,turbine_lat,turb_lon)'相同, –

+1

問題是'lat0'(例如)的類型是'list',因此不支持乘以'float'(雖然乘以int可能不是你想要的!)。在'tunnel_fast'的開頭檢查'print type(lat0)'和'print lat0',你會看到。你以某種方式傳遞錯誤的論點。 –

+0

是的你是對的。在代碼的工作版本中,它是非工作版本中的'和''。 –