def doeqdi(x, y, UP=False):
    """
    Takes digitized x,y, data and returns the dec,inc, assuming an
    equal area projection
    Parameters
    __________________
        x : array of digitized x from point on equal area projection
        y : array of  igitized y from point on equal area projection
        UP : if True, is an upper hemisphere projection
    Output :
        dec : declination
        inc : inclination
    """
    xp, yp = y, x  # need to switch into geographic convention
    r = np.sqrt(xp**2+yp**2)
    z = 1.-r**2
    t = np.arcsin(z)
    if UP == 1:
        t = -t
    p = np.arctan2(yp, xp)
    dec, inc = np.degrees(p) % 360, np.degrees(t)
    return dec, inc