Python numpy.meshgrid() Examples

def generate_anchors_pre(height, width, feat_stride, anchor_scales=(8,16,32), anchor_ratios=(0.5,1,2)):
  """ A wrapper function to generate anchors given different scales
    Also return the number of anchors in variable 'length'
  """
  anchors = generate_anchors(ratios=np.array(anchor_ratios), scales=np.array(anchor_scales))
  A = anchors.shape[0]
  shift_x = np.arange(0, width) * feat_stride
  shift_y = np.arange(0, height) * feat_stride
  shift_x, shift_y = np.meshgrid(shift_x, shift_y)
  shifts = np.vstack((shift_x.ravel(), shift_y.ravel(), shift_x.ravel(), shift_y.ravel())).transpose()
  K = shifts.shape[0]
  # width changes faster, so here it is H, W, C
  anchors = anchors.reshape((1, A, 4)) + shifts.reshape((1, K, 4)).transpose((1, 0, 2))
  anchors = anchors.reshape((K * A, 4)).astype(np.float32, copy=False)
  length = np.int32(anchors.shape[0])

  return anchors, length 

def image_reslice(image, spec, method=None, fill=0, dtype=None, weights=None, image_type=None):
    '''
    image_reslice(image, spec) yields a duplicate of the given image resliced to have the voxels
      indicated by the given image spec. Note that spec may be an image itself.

    Optional arguments that can be passed to image_interpolate() (asside from affine) are allowed
    here and are passed through.
    '''
    if image_type is None and is_image(image): image_type = to_image_type(image)
    spec = to_image_spec(spec)
    image = to_image(image)
    # we make a big mesh and interpolate at these points...
    imsh = spec['image_shape']
    (args, kw) = ([np.arange(n) for n in imsh[:3]], {'indexing': 'ij'})
    ijk = np.asarray([u.flatten() for u in np.meshgrid(*args, **kw)])
    ijk = np.dot(spec['affine'], np.vstack([ijk, np.ones([1,ijk.shape[1]])]))[:3]
    # interpolate here...
    u = image_interpolate(image, ijk, method=method, fill=fill, dtype=dtype, weights=weights)
    return to_image((np.reshape(u, imsh), spec), image_type=image_type) 

def get_point_cloud_from_z(Y, camera_matrix):
  """Projects the depth image Y into a 3D point cloud.
  Inputs:
    Y is ...xHxW
    camera_matrix
  Outputs:
    X is positive going right
    Y is positive into the image
    Z is positive up in the image
    XYZ is ...xHxWx3
  """
  x, z = np.meshgrid(np.arange(Y.shape[-1]),
                     np.arange(Y.shape[-2]-1, -1, -1))
  for i in range(Y.ndim-2):
    x = np.expand_dims(x, axis=0)
    z = np.expand_dims(z, axis=0)
  X = (x-camera_matrix.xc) * Y / camera_matrix.f
  Z = (z-camera_matrix.zc) * Y / camera_matrix.f
  XYZ = np.concatenate((X[...,np.newaxis], Y[...,np.newaxis],
                        Z[...,np.newaxis]), axis=X.ndim)
  return XYZ 

def test_2x3(self):
        # Loading the water depth map
        dat = loadtxt('data/WaterDepth1.dat')
        X, Y = meshgrid(linspace(0., 1000., 50), linspace(0., 1000., 50))
        depth = array(zip(X.flatten(), Y.flatten(), dat.flatten()))
        borders = array([[200, 200], [150, 500], [200, 800], [600, 900], [700, 700], [900, 500], [800, 200], [500, 100], [200, 200]])
        baseline = array([[587.5, 223.07692308], [525., 346.15384615], [837.5, 530.76923077], [525., 530.76923077], [525., 838.46153846], [837.5, 469.23076923]])

        wt_desc = WTDescFromWTG('data/V80-2MW-offshore.wtg').wt_desc
        wt_layout = GenericWindFarmTurbineLayout([WTPC(wt_desc=wt_desc, position=pos) for pos in baseline])

        t = Topfarm(
            baseline_layout = wt_layout,
            borders = borders,
            depth_map = depth,
            dist_WT_D = 5.0,
            distribution='spiral',
            wind_speeds=[4., 8., 20.],
            wind_directions=linspace(0., 360., 36)[:-1]
        )

        t.run()

        self.fail('make save function')
        t.save() 

def generate_anchors(self, image_width: int, image_height: int, num_x_anchors: int, num_y_anchors: int) -> Tensor:
        center_ys = np.linspace(start=0, stop=image_height, num=num_y_anchors + 2)[1:-1]
        center_xs = np.linspace(start=0, stop=image_width, num=num_x_anchors + 2)[1:-1]
        ratios = np.array(self._anchor_ratios)
        ratios = ratios[:, 0] / ratios[:, 1]
        sizes = np.array(self._anchor_sizes)

        # NOTE: it's important to let `center_ys` be the major index (i.e., move horizontally and then vertically) for consistency with 2D convolution
        # giving the string 'ij' returns a meshgrid with matrix indexing, i.e., with shape (#center_ys, #center_xs, #ratios)
        center_ys, center_xs, ratios, sizes = np.meshgrid(center_ys, center_xs, ratios, sizes, indexing='ij')

        center_ys = center_ys.reshape(-1)
        center_xs = center_xs.reshape(-1)
        ratios = ratios.reshape(-1)
        sizes = sizes.reshape(-1)

        widths = sizes * np.sqrt(1 / ratios)
        heights = sizes * np.sqrt(ratios)

        center_based_anchor_bboxes = np.stack((center_xs, center_ys, widths, heights), axis=1)
        center_based_anchor_bboxes = torch.from_numpy(center_based_anchor_bboxes).float()
        anchor_bboxes = BBox.from_center_base(center_based_anchor_bboxes)

        return anchor_bboxes 

def test_dist_sma_radius(self):
        """
        Test that sma and radius values outside of the range have zero probability
        """
        
        for mod in self.allmods:
            if 'dist_sma_radius' in mod.__dict__:
                with RedirectStreams(stdout=self.dev_null):
                    pp = mod(**self.spec)
                
                a = np.logspace(np.log10(pp.arange[0].value/10.),np.log10(pp.arange[1].value*100),100)
                Rp = np.logspace(np.log10(pp.Rprange[0].value/10.),np.log10(pp.Rprange[1].value*100),100)
                
                aa, RR = np.meshgrid(a,Rp)
                
                fr = pp.dist_sma_radius(aa,RR)
                self.assertTrue(np.all(fr[aa < pp.arange[0].value] == 0),'dist_sma_radius low bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[aa > pp.arange[1].value] == 0),'dist_sma_radius high bound failed on sma for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR < pp.Rprange[0].value] == 0),'dist_sma_radius low bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[RR > pp.Rprange[1].value] == 0),'dist_sma_radius high bound failed on radius for %s'%mod.__name__)
                self.assertTrue(np.all(fr[(aa > pp.arange[0].value) & (aa < pp.arange[1].value) & (RR > pp.Rprange[0].value) & (RR < pp.Rprange[1].value)] > 0),'dist_sma_radius is improper pdf for %s'%mod.__name__) 

def calcfbetaInput(self):
        # table 17 in Leinert et al. (1998)
        # Zodiacal Light brightness function of solar LON (rows) and LAT (columns)
        # values given in W m−2 sr−1 μm−1 for a wavelength of 500 nm
        path = os.path.split(inspect.getfile(self.__class__))[0]
        Izod = np.loadtxt(os.path.join(path, 'Leinert98_table17.txt'))*1e-8 # W/m2/sr/um
        # create data point coordinates
        lon_pts = np.array([0., 5, 10, 15, 20, 25, 30, 35, 40, 45, 60, 75, 90,
                105, 120, 135, 150, 165, 180]) # deg
        lat_pts = np.array([0., 5, 10, 15, 20, 25, 30, 45, 60, 75, 90]) # deg
        y_pts, x_pts = np.meshgrid(lat_pts, lon_pts)
        points = np.array(list(zip(np.concatenate(x_pts), np.concatenate(y_pts))))
        # create data values, normalized by (90,0) value
        z = Izod/Izod[12,0]
        values = z.reshape(z.size)
        return  points, values 

def get_L_cluster_cut(L, node_label):
  adj = L - np.diag(np.diag(L))
  adj[adj != 0] = 1.0
  num_nodes = adj.shape[0]
  idx_row, idx_col = np.meshgrid(range(num_nodes), range(num_nodes))
  idx_row, idx_col = idx_row.flatten().astype(
      np.int64), idx_col.flatten().astype(np.int64)
  mask = (node_label[idx_row] == node_label[idx_col]).reshape(
      num_nodes, num_nodes).astype(np.float)

  adj_cluster = adj * mask
  adj_cut = adj - adj_cluster
  L_cut = get_laplacian(adj_cut, graph_laplacian_type='L4')
  L_cluster = get_laplacian(adj_cluster, graph_laplacian_type='L4')

  return L_cluster, L_cut 

def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 

def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 

def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 

def ov_order(model, slc=None):
    nocc = k_nocc(model)
    if slc is None:
        slc = numpy.ones((len(model.mo_coeff), model.mo_coeff[0].shape[1]), dtype=bool)
    e_occ = tuple(e[:o][s[:o]] for e, o, s in zip(model.mo_energy, nocc, slc))
    e_virt = tuple(e[o:][s[o:]] for e, o, s in zip(model.mo_energy, nocc, slc))
    sort_o = []
    sort_v = []
    for o in e_occ:
        for v in e_virt:
            _v, _o = numpy.meshgrid(v, o)
            sort_o.append(_o.reshape(-1))
            sort_v.append(_v.reshape(-1))
    sort_o, sort_v = numpy.concatenate(sort_o), numpy.concatenate(sort_v)
    vals = numpy.array(
        list(zip(sort_o, sort_v)),
        dtype=[('o', sort_o[0].dtype), ('v', sort_v[0].dtype)]
    )
    result = numpy.argsort(vals, order=('o', 'v'))
    # Double for other blocks
    return numpy.concatenate([result, result + len(result)]) 

def test_draw_polygon_2D(self):
        X, Y = np.meshgrid(np.linspace(-8, 8, 25, dtype=float), np.linspace(-8, 8, 25, dtype=float))
        nodes = np.vstack([X.reshape(-1), Y.reshape(-1)]).T
        tri = scipy.spatial.Delaunay(nodes)
        poly = np.array([[5.5, 5.5], [5.5, -5.5], [-5.5, -5.5], [-5.5, 5.5]])
        #angles = np.linspace(0, 2*np.pi, endpoint=False, num=20)
        #poly = np.vstack([5*np.sin(angles), 5*np.cos(angles)]).T
        trs = tri.simplices[:, [1,0,2]]
        new_points, _ = electrode_placement._draw_polygon_2D(
            poly, tri.points, trs, ends=True)
        bar = np.mean(new_points[trs], axis=1)
        m = new_points[trs[:, 1:]] -\
            new_points[trs[:, 0]][:, None, :]
        area = .5 * -np.linalg.det(m)
        inside = electrode_placement._point_inside_polygon(poly, bar, tol=1e-3)
        #plt.triplot(new_points[:, 0], new_points[:, 1], tri.simplices.copy())
        #plt.show()
        assert np.isclose(np.sum(area[inside]), 11**2) 

def test_inside_complex_polygon(self):
        X, Y = np.meshgrid(np.linspace(-8, 8, 17, dtype=float), np.linspace(-8, 8, 17, dtype=float))
        nodes = np.vstack([X.reshape(-1), Y.reshape(-1)]).T
        tri = scipy.spatial.Delaunay(nodes)
        poly = np.array([[5, 5], [5, -5], [-5, -5], [-5, 5]])
        hole1 = np.array([[2, 2], [2, -2], [-2, -2], [-2, 2]])
        hole2 = np.array([[4, 4], [4, 3], [3, 3], [3, 4]])
        trs = tri.simplices[:, [1,0,2]]
        inside = electrode_placement._inside_complex_polygon(poly, tri.points,
                                                             trs,
                                                             holes=[hole1, hole2],
                                                             tol=1e-3)
        m = tri.points[trs[inside, 1:]] -\
            tri.points[trs[inside, 0]][:, None, :]
        area = .5 * -np.linalg.det(m)
        assert np.isclose(np.sum(area), 100-16-1) 

def test_find_tetrahedron_with_points_non_convex(self, sphere3_msh):
        X, Y, Z = np.meshgrid(np.linspace(-100, 100, 100),
                              np.linspace(-40, 40, 100), [0])
        points = np.vstack([X.reshape(-1), Y.reshape(-1), Z.reshape(-1)]).T
        dist = np.linalg.norm(points, axis=1)

        msh = sphere3_msh.crop_mesh(5)
        points_outside = points[(dist > 95) + (dist < 89)]
        th_with_points, bar = msh.find_tetrahedron_with_points(points_outside)
        assert np.all(th_with_points == -1)
        assert np.allclose(bar, 0)

        points_inside = points[(dist <= 94) * (dist >= 91)]
        th_with_points, bar = msh.find_tetrahedron_with_points(points_inside)
        eps = 1e-3
        assert np.all(th_with_points != -1)
        assert np.all(bar >= 0 - eps)
        assert np.all(bar <= 1. + eps)
        th_coords = \
            msh.nodes[msh.elm[th_with_points]]
        assert np.allclose(np.einsum('ikj, ik -> ij', th_coords, bar), points_inside) 

def test_inside_volume(self, sphere3_msh):
        X, Y, Z = np.meshgrid(np.linspace(-100, 100, 100),
                              np.linspace(-40, 40, 10), [0])
        np.random.seed(0)
        points = np.vstack([X.reshape(-1), Y.reshape(-1), Z.reshape(-1)]).T
        points += np.random.random_sample(points.shape) - .5
        dist = np.linalg.norm(points, axis=1)

        msh = sphere3_msh.crop_mesh([4, 5])
        points_outside = points[(dist > 95) + (dist < 84)]
        inside = msh.test_inside_volume(points_outside)
        assert np.all(~inside)

        points_inside = points[(dist <= 94) * (dist >= 86)]
        inside = msh.test_inside_volume(points_inside)
        assert np.all(inside) 

def test_to_nonlinear_grid(self, sphere3_msh):
        import nibabel
        data = sphere3_msh.elm.tag1
        f = mesh_io.ElementData(data, mesh=sphere3_msh)
        affine = np.array([[1, 0, 0, -100.5],
                           [0, 1, 0, -5],
                           [0, 0, 1, 0],
                           [0, 0, 0, 1]], dtype=float)
        x, y, z = np.meshgrid(np.arange(-100, 100),
                              np.arange(-5, 5),
                              np.arange(0, 1),
                              indexing='ij')
        nonl_transform = np.concatenate(
            (x[..., None], y[..., None], z[..., None]), axis=3).astype(float)
        img = nibabel.Nifti1Pair(nonl_transform, affine)
        tempf = tempfile.NamedTemporaryFile(suffix='.nii', delete=False)
        fn = tempf.name
        tempf.close()
        nibabel.save(img, fn)
        interp = f.to_deformed_grid(fn, fn, method='assign')
        interp = interp.get_data()
        assert np.isclose(interp[100, 5, 0], 3)
        assert np.isclose(interp[187, 5, 0], 4)
        assert np.isclose(interp[193, 5, 0], 5)
        assert np.isclose(interp[198, 5, 0], 0) 

def test_to_nonlinear_grid_crop(self, sphere3_msh):
        import nibabel
        data = sphere3_msh.elm.tag1
        f = mesh_io.ElementData(data, mesh=sphere3_msh)
        affine = np.array([[1, 0, 0, -100.5],
                           [0, 1, 0, -5],
                           [0, 0, 1, 0],
                           [0, 0, 0, 1]], dtype=float)
        x, y, z = np.meshgrid(np.arange(-100, 100),
                              np.arange(-5, 5),
                              np.arange(0, 1),
                              indexing='ij')
        nonl_transform = np.concatenate(
            (x[..., None], y[..., None], z[..., None]), axis=3).astype(float)
        img = nibabel.Nifti1Pair(nonl_transform, affine)
        tempf = tempfile.NamedTemporaryFile(suffix='.nii', delete=False)
        fn = tempf.name
        tempf.close()
        nibabel.save(img, fn)
        interp = f.to_deformed_grid(fn, fn, tags=3, method='assign')
        interp = interp.get_data()
        assert np.isclose(interp[100, 5, 0], 3)
        assert np.isclose(interp[187, 5, 0], 0)
        assert np.isclose(interp[193, 5, 0], 0)
        assert np.isclose(interp[198, 5, 0], 0) 

def test_from_data_grid_vec_extra_args(self, sphere3_msh):
        X, Y, Z = np.meshgrid(np.linspace(-100, 100, 201),
                              np.linspace(-100, 100, 201),
                              np.linspace(-1, 1, 3),
                              indexing='ij')
        V = np.stack([X, Y, Z], axis=3)
        affine = np.array([[1, 0, 0, -100],
                           [0, 1, 0, -100],
                           [0, 0, 1, -1],
                           [0, 0, 0, 1]], dtype=float)
        ed = mesh_io.ElementData.from_data_grid(sphere3_msh, V, affine,
                                             cval=0.0, order=1)
        bar = sphere3_msh.elements_baricenters().value
        in_area = (bar[:, 2] >= -1) * (bar[:, 2] <= 1)
        assert np.allclose(bar[in_area], ed.value[in_area])
        assert np.allclose(ed.value[~in_area], 0) 

def test_big_cell():
    import time

    a = 1
    ncell = (2, 2, 2)
    Lvecs = np.diag(ncell) * a
    unit_cell = np.zeros((4, 3))
    unit_cell[1:] = (np.ones((3, 3)) - np.eye(3)) * a / 2

    grid = np.meshgrid(*map(np.arange, ncell), indexing="ij")
    shifts = np.stack(list(map(np.ravel, grid)), axis=1)
    supercell = (shifts[:, np.newaxis] + unit_cell[np.newaxis]).reshape(1, -1, 3)

    configs = supercell.repeat(1000, axis=0)
    configs += np.random.randn(*configs.shape) * 0.1

    df = run(Lvecs, configs, 8)
    df = df.groupby("qmag").mean().reset_index()

    large_q = df[-35:-10]["Sq"]
    mean = np.mean(large_q - 1)
    rms = np.sqrt(np.mean((large_q - 1) ** 2))
    assert np.abs(mean) < 0.01, mean
    assert rms < 0.1, rms 

def plotc(parameters, ax=[], color='k', gridsize=(101, 101)):
    """
    Plot a linear classifier in a 2D scatterplot.

    INPUT   (1) tuple 'parameters': consists of a list of class proportions
                (1 by K classes), an array of class means (K classes by
                D features), an array of class-covariance matrices (D features
                by D features by K classes)
            (2) object 'ax': axes of a pyplot figure or subject (def: empty)
            (3) str 'colors': colors of the contours in the plot (def: 'k')
            (4) tuple 'gridsize': number of points in the grid
                (def: (101, 101))
    OUTPUT  None
    """
    # Check for figure object
    if fig:
        ax = fig.gca()
    else:
        fig, ax = plt.subplots()

    # Get axes limits
    xl = ax.get_xlim()
    yl = ax.get_ylim()

    # Define grid
    gx = np.linspace(xl[0], xl[1], gridsize[0])
    gy = np.linspace(yl[0], yl[1], gridsize[1])
    x, y = np.meshgrid(gx, gy)
    xy = np.vstack((x.ravel(), y.ravel())).T

    # Values of grid
    z = np.dot(xy, parameters[:-1, :]) + parameters[-1, :]
    z = np.reshape(z[:, 0] - z[:, 1], gridsize)

    # Plot grid
    ax.contour(x, y, z, levels=0, colors=colors) 

def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up BH query object
    print('Loading BH map...')
    bh = BHQuery()

    print('Querying map...')
    ebv = bh.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'bh.png'
    img.save(fname)

    return 0 

def main():
    w,h = (2*2056, 2*int(2056*(20./200.)))
    l_0 = 0.

    # Set up MarshallQuery object
    print('Loading Marshall map...')
    query = MarshallQuery()

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-100.+l_0, 100.+l_0, 2*w)
    b = np.linspace(-10., 10., 2*h)
    dl = l[1] - l[0]
    db = b[1] - b[0]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * dl
    b += (np.random.random(l.shape) - 0.5) * db

    A = np.empty(l.shape+(3,), dtype='f8')

    for k,d in enumerate([1., 2.5, 5.]):
        coords = SkyCoord(l*u.deg, b*u.deg, d*u.kpc, frame='galactic')

        # Get the mean dust extinction at each coordinate
        print('Querying map...')
        A[:,:,k] = query(coords, return_sigma=False)

    A[:,:,2] -= A[:,:,1]
    A[:,:,1] -= A[:,:,0]

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(A[::-1,::-1,:], 0., 1., fill=255)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'marshall.png'
    img.save(fname)

    return 0 

def main():
    w,h = (2056,1024)
    l_0 = 130.

    # Set up Bayestar query object
    print('Loading bayestar map...')
    bayestar = BayestarQuery(max_samples=1)

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    ebv = np.empty(l.shape+(3,), dtype='f8')

    for k,d in enumerate([0.5, 1.5, 5.]):
        # d = 5.    # We'll query integrated reddening to a distance of 5 kpc
        coords = SkyCoord(l*u.deg, b*u.deg, d*u.kpc, frame='galactic')

        # Get the dust median reddening at each coordinate
        print('Querying map...')
        ebv[:,:,k] = bayestar.query(coords, mode='median')

    ebv[:,:,2] -= ebv[:,:,1]
    ebv[:,:,1] -= ebv[:,:,0]

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1,:], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'bayestar.png'
    img.save(fname)

    return 0 

def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up Lenz+(2017) query object
    print('Loading Lenz+(2017) map...')
    q = Lenz2017Query()

    print('Querying map...')
    ebv = q.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 0.05)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'lenz2017.png'
    img.save(fname)

    return 0 

def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    planck_components = [
        ('ebv', 0., 1.5),
        ('radiance', 0., 1.5),
        ('tau', 0., 1.5),
        ('temp', 15.*u.K, 25.*u.K),
        ('err_temp', 0.*u.K, 4.*u.K),
        ('beta', 1., 3.),
        ('err_beta', 0., 0.2)]

    for component,vmin,vmax in planck_components:
        # Set up Planck query object
        print('Loading Planck map...')
        planck = PlanckQuery(component=component)

        print('Querying map...')
        res = planck.query(coords)

        # Convert the output array to a PIL image and save
        print('Saving image...')
        img = numpy2pil(res[::-1,::-1], vmin, vmax)
        img = img.resize((w,h), resample=PIL.Image.LANCZOS)
        fname = 'planck_{}.png'.format(component)
        img.save(fname)

    return 0 

def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up SFD query object
    print('Loading SFD map...')
    sfd = SFDQuery()

    print('Querying map...')
    ebv = sfd.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'sfd.png'
    img.save(fname)

    return 0 

def main():
    w,h = (2*2056, 2*int(2056*(30./200.)))
    l_0 = 122.5

    # Set up IPHASquery object
    print('Loading IPHAS map...')
    iphas = IPHASQuery()

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-100.+l_0, 100.+l_0, 2*w)
    b = np.linspace(-15., 15., 2*h)
    dl = l[1] - l[0]
    db = b[1] - b[0]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * dl
    b += (np.random.random(l.shape) - 0.5) * db

    A = np.empty(l.shape+(3,), dtype='f8')

    for k,d in enumerate([0.5, 1.5, 5.]):
        # d = 5.    # We'll query integrated reddening to a distance of 5 kpc
        coords = SkyCoord(l*u.deg, b*u.deg, d*u.kpc, frame='galactic')

        # Get the dust median reddening at each coordinate
        print('Querying map...')
        A[:,:,k] = iphas.query(coords, mode='random_sample')

    A[:,:,2] -= A[:,:,1]
    A[:,:,1] -= A[:,:,0]

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(A[::-1,::-1,:], 0., 4.5, fill=255)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'iphas.png'
    img.save(fname)

    return 0 

def show_trace_2d(f, results):
    plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = np.meshgrid(np.arange(-5.5, 1.0, 0.1), np.arange(-3.0, 1.0, 0.1))
    plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    plt.xlabel('x1')
    plt.ylabel('x2')
    plt.show() 

def show_trace_2d(f, results):
    plt.plot(*zip(*results), '-o', color='#ff7f0e')
    x1, x2 = np.meshgrid(np.arange(-5.5, 1.0, 0.1), np.arange(-3.0, 1.0, 0.1))
    plt.contour(x1, x2, f(x1, x2), colors='#1f77b4')
    plt.xlabel('x1')
    plt.ylabel('x2')
    plt.show()