pbo_util.py 文件源码

def getROIstations(geo_point,radiusParam,data,header):
    '''
    This function returns the 4ID station codes for the stations in a region

    The region of interest is defined by the geographic coordinate and a window size

    @param geo_point: The geographic (lat,lon) coordinate of interest
    @param radiusParam: An overloaded radius of interest [km] or latitude and longitude window [deg] around the geo_point
    @param data: Stabilized (or unstabilized) data generated from the data fetcher or out of stab_sys
    @param header: Header dictionary with stations metadata keyed by their 4ID code. This is output with the data.

    @return station_list, list of site 4ID codes in the specified geographic region
     '''
    ccPos = (geo_point[0]*np.pi/180, geo_point[1]*np.pi/180)
    if np.isscalar(radiusParam):
        station_list = []
        for ii in header.keys():
            coord = (header[ii]['refNEU'][0]*np.pi/180,(header[ii]['refNEU'][1]-360)*np.pi/180)
            dist = 6371*2*np.arcsin(np.sqrt(np.sin((ccPos[0]-coord[0])/2)**2+np.cos(ccPos[0])*np.cos(coord[0])*np.sin((ccPos[1]-coord[1])/2)**2))
            if np.abs(dist) < radiusParam:
                station_list.append(header[ii]['4ID'])
    else:
        # overloaded radiusParam term to be radius or lat/lon window size
        latWin = radiusParam[0]/2
        lonWin = radiusParam[1]/2
        station_list = []

        try:
            for ii in header.keys():
                coord = (header[ii]['refNEU'][0],(header[ii]['refNEU'][1]-360))
                if (geo_point[0]-latWin)<=coord[0]<=(geo_point[0]+latWin) and (geo_point[1]-lonWin)<=coord[1]<=(geo_point[1]+lonWin):
                    station_list.append(header[ii]['4ID'])
        except:
            station_list = None

    return station_list