Get lons and lats from x, y in projection coordinates.
972159
def _lonlat_from_geos_angle(x, y, geos_area):
    """"""
    h = (geos_area.proj_dict['h'] + geos_area.proj_dict['a']) / 1000
    b__ = (geos_area.proj_dict['a'] / geos_area.proj_dict['b']) ** 2

    sd = np.sqrt((h * np.cos(x) * np.cos(y)) ** 2 -
                 (np.cos(y)**2 + b__ * np.sin(y)**2) *
                 (h**2 - (geos_area.proj_dict['a'] / 1000)**2))
    # sd = 0

    sn = (h * np.cos(x) * np.cos(y) - sd) / (np.cos(y)**2 + b__ * np.sin(y)**2)
    s1 = h - sn * np.cos(x) * np.cos(y)
    s2 = sn * np.sin(x) * np.cos(y)
    s3 = -sn * np.sin(y)
    sxy = np.sqrt(s1**2 + s2**2)

    lons = np.rad2deg(np.arctan2(s2, s1)) + geos_area.proj_dict.get('lon_0', 0)
    lats = np.rad2deg(-np.arctan2(b__ * s3, sxy))

    return lons, lats
286849
def get_lons_from_cartesian(x__, y__):
    
    return rad2deg(arccos(x__ / sqrt(x__ ** 2 + y__ ** 2))) * sign(y__)
231558
def xyz2lonlat(x, y, z):
    """"""
    lon = xu.rad2deg(xu.arctan2(y, x))
    lat = xu.rad2deg(xu.arctan2(z, xu.sqrt(x**2 + y**2)))
    return lon, lat
572333
def lon_lat(value):
    
    lon, lat = value.split()
    return longitude(lon), latitude(lat)
352048
def native_projection_explanation(ax):
    """"""
    ax.set_title('Longitude and Latitude', size=18, y=1.1)

    # Hide the azimuth labels
    ax.set_azimuth_ticklabels([])

    # Make the axis tick labels visible:
    ax.set_xticks(np.radians(np.arange(-80, 90, 10)))
    ax.tick_params(label1On=True)

    ax.set_xlabel('Longitude')

    xlabel_halo(ax)
    return ax
630893
def interp_value(self, lat, lon, indexed=False):
        
        (px, py) = self.grid_coordinates.projection_to_raster_coords(lat, lon)
        if indexed:
            return self.raster_data[round(py), round(px)]
        else:
#             from scipy.interpolate import interp2d
#             f_interp = interp2d(self.grid_coordinates.x_axis, self.grid_coordinates.y_axis, self.raster_data, bounds_error=True)
#             return f_interp(lon, lat)[0]
            from scipy.ndimage import map_coordinates
            ret = map_coordinates(self.raster_data, [[py], [px]], order=1)  # linear interp
            return ret[0]
533893
def from_radians(cls, lat_radians, long_radians):
        '''
        
        '''
        return cls(math.degrees(lat_radians), math.degrees(long_radians))
511776
def raster_to_projection_coords(self, pixel_x, pixel_y):
        
        h_px_py = np.array([1, pixel_x, pixel_y])
        gt = np.array([[1, 0, 0], self.geotransform[0:3], self.geotransform[3:6]])
        arr = np.inner(gt, h_px_py)
        return arr[2], arr[1]
1281742
def get_xyz_1d(self, xcoord, x, ycoord, y, u, v):
        """"""
        xclose = xcoord.indexes[xcoord.name].get_loc(x, method='nearest')
        yclose = ycoord.indexes[ycoord.name].get_loc(y, method='nearest')
        uval = u[yclose, xclose].values
        vval = v[yclose, xclose].values
        return xcoord[xclose].values, ycoord[yclose].values, uval, vval
903802
def get_pixel_indices(self, lats, lons):
        
        theta = np.radians(90. - lats)
        phi = np.radians(lons)
        return hp.ang2pix(self.nside, theta, phi, self.nest)
245225
def area2lonlat(dataarray):
    """"""
    area = dataarray.attrs['area']
    lons, lats = area.get_lonlats_dask()
    lons = xr.DataArray(lons, dims=['y', 'x'],
                        attrs={'name': "longitude",
                               'standard_name': "longitude",
                               'units': 'degrees_east'},
                        name='longitude')
    lats = xr.DataArray(lats, dims=['y', 'x'],
                        attrs={'name': "latitude",
                               'standard_name': "latitude",
                               'units': 'degrees_north'},
                        name='latitude')
    dataarray.attrs['coordinates'] = 'longitude latitude'

    return [dataarray, lons, lats]
1081224
def format_coord(self, x, y):
        """"""
        p, b = stereonet_math.geographic2plunge_bearing(x, y)
        s, d = stereonet_math.geographic2pole(x, y)
        pb = u'P/B={:0.0f}\u00b0/{:03.0f}\u00b0'.format(p[0], b[0])
        sd = u'S/D={:03.0f}\u00b0/{:0.0f}\u00b0'.format(s[0], d[0])
        return u'{}, {}'.format(pb, sd)
880685
def _get_projection(self):
        """"""
        try:
            proj_str = self.nc.attrs['gdal_projection']
        except TypeError:
            proj_str = self.nc.attrs['gdal_projection'].decode()

        # Check the a/b/h units
        radius_a = proj_str.split('+a=')[-1].split()[0]
        if float(radius_a) > 10e3:
            units = 'm'
            scale = 1.0
        else:
            units = 'km'
            scale = 1e3

        if 'units' not in proj_str:
            proj_str = proj_str + ' +units=' + units

        area_extent = (float(self.nc.attrs['gdal_xgeo_up_left']) / scale,
                       float(self.nc.attrs['gdal_ygeo_low_right']) / scale,
                       float(self.nc.attrs['gdal_xgeo_low_right']) / scale,
                       float(self.nc.attrs['gdal_ygeo_up_left']) / scale)

        return proj_str, area_extent
941549
def lat_to_deg(lat):
    """"""
    if isinstance(lat, str) and (':' in lat):
        # TODO: handle other coordinate systems
        lat_deg = dmsStrToDeg(lat)
    else:
        lat_deg = float(lat)
    return lat_deg
525624
def get_full_lonlats(self):
        
        if self.lons is not None and self.lats is not None:
            return self.lons, self.lats

        self.lons, self.lats = self._get_full_lonlats()
        self.lons = da.from_delayed(self.lons, dtype=self["EARTH_LOCATIONS"].dtype,
                                    shape=(self.scanlines, self.pixels))
        self.lats = da.from_delayed(self.lats, dtype=self["EARTH_LOCATIONS"].dtype,
                                    shape=(self.scanlines, self.pixels))
        return self.lons, self.lats
472469
def from_points(lons, lats):
        
        north, west = max(lats), min(lons)
        south, east = min(lats), max(lons)
        return BoundingBox(north=north, west=west, south=south, east=east)
1119480
def _load_GeoTransform(self):
        """"""
        def load_lon():
            return arange(ds.RasterXSize)*b[1]+b[0]

        def load_lat():
            return arange(ds.RasterYSize)*b[5]+b[3]
        ds = self.ds
        b = self.ds.GetGeoTransform()  # bbox, interval
        if with_dask:
            lat = Array(
                {('lat', 0): (load_lat,)}, 'lat', (self.ds.RasterYSize,),
                shape=(self.ds.RasterYSize,), dtype=float)
            lon = Array(
                {('lon', 0): (load_lon,)}, 'lon', (self.ds.RasterXSize,),
                shape=(self.ds.RasterXSize,), dtype=float)
        else:
            lat = load_lat()
            lon = load_lon()
        return Variable(('lat',), lat), Variable(('lon',), lon)