<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def gen_point_query( vectors, raster, band=1, layer=0, nodata=None, affine=None, interpolate='bilinear', property_name='value', geojson_out=False):
 """ Given a set of vector features and a raster, generate raster values at each vertex of the geometry For features with point geometry, the values will be a 1D with the index refering to the feature For features with other geometry types, it effectively creates a 2D list, such that the first index is the feature, the second is the vertex within the geometry Parameters vectors: path to an vector source or geo-like python objects raster: ndarray or path to a GDAL raster source If ndarray is passed, the `transform` kwarg is required. layer: int or string, optional If `vectors` is a path to an fiona source, specify the vector layer to use either by name or number. defaults to 0 band: int, optional If `raster` is a GDAL source, the band number to use (counting from 1). defaults to 1. nodata: float, optional If `raster` is a GDAL source, this value overrides any NODATA value specified in the file's metadata. If `None`, the file's metadata's NODATA value (if any) will be used. defaults to `None`. affine: Affine instance required only for ndarrays, otherwise it is read from src interpolate: string 'bilinear' or 'nearest' interpolation property_name: string name of property key if geojson_out geojson_out: boolean generate GeoJSON-like features (default: False) original feature geometry and properties will be retained point query values appended as additional properties. Returns ------- generator of arrays (if ``geojson_out`` is False) generator of geojson features (if ``geojson_out`` is True) """

    if interpolate not in ['nearest', 'bilinear']:
        raise ValueError("interpolate must be nearest or bilinear")

    features_iter = read_features(vectors, layer)

    with Raster(raster, nodata=nodata, affine=affine, band=band) as rast:

        for feat in features_iter:
            geom = shape(feat['geometry'])
            vals = []
            for x, y in geom_xys(geom):
                if interpolate == 'nearest':
                    r, c = rast.index(x, y)
                    window = ((int(r), int(r+1)), (int(c), int(c+1)))
                    src_array = rast.read(window=window, masked=True).array
                    val = src_array[0, 0]
                    if val is masked:
                        vals.append(None)
                    else:
                        vals.append(asscalar(val))

                elif interpolate == 'bilinear':
                    window, unitxy = point_window_unitxy(x, y, rast.affine)
                    src_array = rast.read(window=window, masked=True).array
                    vals.append(bilinear(src_array, *unitxy))

            if len(vals) == 1:
                vals = vals[0]  # flatten single-element lists

            if geojson_out:
                if 'properties' not in feat:
                    feat['properties'] = {}
                feat['properties'][property_name] = vals
                yield feat
            else:
                yield vals