Python shapely.wkt.loads() Examples

The following are 30 code examples of shapely.wkt.loads(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module shapely.wkt , or try the search function .
Example #1
Source File: test_encoder.py    From mapbox-vector-tile with MIT License 8 votes vote down vote up
def test_invalid_geometry_make_valid(self):
        from mapbox_vector_tile import encode
        from mapbox_vector_tile.encoder import on_invalid_geometry_make_valid
        import shapely.geometry
        import shapely.wkt
        geometry = 'POLYGON ((10 10, 20 10, 20 20, 15 15, 15 5, 10 10))'
        shape = shapely.wkt.loads(geometry)
        self.assertFalse(shape.is_valid)
        feature = dict(geometry=shape, properties={})
        source = dict(name='layername', features=[feature])
        pbf = encode(source,
                     on_invalid_geometry=on_invalid_geometry_make_valid)
        result = decode(pbf)
        self.assertEqual(1, len(result['layername']['features']))
        valid_geometry = result['layername']['features'][0]['geometry']
        self.assertTrue(valid_geometry['type'], 'MultiPolygon')
        multipolygon = shapely.geometry.shape(valid_geometry)
        self.assertTrue(multipolygon.geom_type, 'MultiPolygon')
        self.assertTrue(multipolygon.is_valid)

        for poly in multipolygon.geoms:
            self.assertTrue(poly.is_valid) 
Example #2
Source File: elements.py    From momepy with MIT License 7 votes vote down vote up
def _densify(self, geom, segment):
        """
        Returns densified geoemtry with segments no longer than `segment`.
        """
        # temporary solution for readthedocs fail. - cannot mock osgeo
        try:
            from osgeo import ogr
        except ModuleNotFoundError:
            import warnings

            warnings.warn("OGR (GDAL) is required.")

        poly = geom
        wkt = geom.wkt  # shapely Polygon to wkt
        geom = ogr.CreateGeometryFromWkt(wkt)  # create ogr geometry
        geom.Segmentize(segment)  # densify geometry by 2 metres
        geom.CloseRings()  # fix for GDAL 2.4.1 bug
        wkt2 = geom.ExportToWkt()  # ogr geometry to wkt
        try:
            new = loads(wkt2)  # wkt to shapely Polygon
            return new
        except Exception:
            return poly 
Example #3
Source File: meta.py    From gbdxtools with MIT License 6 votes vote down vote up
def _parse_geoms(self, **kwargs):
        """ Finds supported geometry types, parses them and returns the bbox """
        bbox = kwargs.get('bbox', None)
        wkt_geom = kwargs.get('wkt', None)
        geojson = kwargs.get('geojson', None)
        if bbox is not None:
            g = box(*bbox)
        elif wkt_geom is not None:
            g = wkt.loads(wkt_geom)
        elif geojson is not None:
            g = shape(geojson)
        else:
            return None
        if self.proj is None:
            return g
        else:
            return self._reproject(g, from_proj=kwargs.get('from_proj', 'EPSG:4326')) 
Example #4
Source File: bench_encode.py    From mapbox-vector-tile with MIT License 6 votes vote down vote up
def make_layers(shapes, geom_dicts=False):
    print("Creating layers with 10 shapes each")
    layers = []
    i = 0
    features = []
    for shape in shapes:
        try:
            geom = loads_wkt(shape.strip())
            if geom_dicts:
                feature = {"geometry": mapping(geom), "properties": {}}
            else:
                feature = {"geometry": geom, "properties": {}}
            features.append(feature)
            if i >= 10:
                layers.append(features)
                features = []
                i = 0
            i += 1
        except:
            pass
    layers.append(features)
    return layers 
Example #5
Source File: test_encoder.py    From mapbox-vector-tile with MIT License 6 votes vote down vote up
def test_bowtie_self_touching(self):
        from mapbox_vector_tile import encode
        from mapbox_vector_tile.encoder import on_invalid_geometry_make_valid
        import shapely.geometry
        import shapely.wkt
        bowtie = ('POLYGON ((0 0, 0 2, 1 1, 2 2, 2 0, 1 1, 0 0))')
        shape = shapely.wkt.loads(bowtie)
        self.assertFalse(shape.is_valid)
        feature = dict(geometry=shape, properties={})
        source = dict(name='layername', features=[feature])
        pbf = encode(source,
                     on_invalid_geometry=on_invalid_geometry_make_valid)
        result = decode(pbf)
        self.assertEqual(1, len(result['layername']['features']))
        valid_geometries = result['layername']['features'][0]['geometry']
        multipolygon = shapely.geometry.shape(valid_geometries)
        self.assertEqual(2, len(multipolygon.geoms))
        shape1, shape2 = multipolygon.geoms
        self.assertTrue(shape1.is_valid)
        self.assertTrue(shape2.is_valid)
        self.assertGreater(shape1.area, 0)
        self.assertGreater(shape2.area, 0) 
Example #6
Source File: test_indexed_geometries.py    From maup with MIT License 6 votes vote down vote up
def test_intersections_correct_when_all_overlapping(four_square_grid, square):
    indexed = IndexedGeometries(four_square_grid)
    overlaps = indexed.intersections(square)

    expected_polygons = [
        wkt.loads(p)
        for p in [
            "POLYGON ((0.5 1, 1 1, 1 0.5, 0.5 0.5, 0.5 1))",
            "POLYGON ((1 1.5, 1 1, 0.5 1, 0.5 1.5, 1 1.5))",
            "POLYGON ((1 0.5, 1 1, 1.5 1, 1.5 0.5, 1 0.5))",
            "POLYGON ((1 1, 1 1.5, 1.5 1.5, 1.5 1, 1 1))",
        ]
    ]

    for p in expected_polygons:
        assert any(overlap == p for overlap in overlaps)

    for p in overlaps:
        assert any(p == expected for expected in expected_polygons) 
Example #7
Source File: test_encoder.py    From mapbox-vector-tile with MIT License 6 votes vote down vote up
def test_validate_generates_rounding_error(self):
        from mapbox_vector_tile import encode
        from mapbox_vector_tile.encoder import on_invalid_geometry_make_valid
        import shapely.geometry
        import shapely.wkt
        bowtie = ('POLYGON((0 0, 1 1, 0 1, 1 0, 0 0))')
        shape = shapely.wkt.loads(bowtie)
        self.assertFalse(shape.is_valid)
        feature = dict(geometry=shape, properties={})
        source = dict(name='layername', features=[feature])
        pbf = encode(source,
                     on_invalid_geometry=on_invalid_geometry_make_valid)
        result = decode(pbf)
        features = result['layername']['features']
        self.assertEqual(1, len(features))
        self.assertEqual(features[0]['geometry']['type'], 'Polygon')
        shape = shapely.geometry.shape(features[0]['geometry'])
        self.assertEqual(shape.geom_type, 'Polygon')
        self.assertTrue(shape.is_valid)
        self.assertGreater(shape.area, 0) 
Example #8
Source File: test_encoder.py    From mapbox-vector-tile with MIT License 6 votes vote down vote up
def test_flipped_geometry_produces_multipolygon(self):
        from mapbox_vector_tile import encode
        from mapbox_vector_tile.encoder import on_invalid_geometry_make_valid
        import shapely.wkt
        shape = shapely.wkt.loads('POLYGON ((3449 1939, 3476 1967, 3473 1996, 3483 2027, 3542 2119, 3538 2160, 3563 2233, 3602 2255, 3639 2326, 3629 2388, 3573 2455, 3594 2493, 3558 2533, 3573 2549, 3518 2572, 3502 2592, 3505 2607, 3513 2614, 3535 2616, 3537 2610, 3535 2602, 3537 2599, 3548 2607, 3551 2636, 3528 2634, 3537 2668, 3549 2670, 3528 2711, 3550 2667, 3532 2635, 3550 2641, 3553 2613, 3549 2602, 3540 2596, 3512 2610, 3506 2589, 3576 2552, 3576 2543, 3563 2535, 3596 2506, 3597 2494, 3587 2469, 3589 2451, 3636 2385, 3644 2326, 3605 2251, 3566 2230, 3547 2122, 3482 2014, 3479 1966, 3455 1944, 3458 1910, 3449 1902, 3449 1939))')  # noqa
        features = [dict(geometry=shape, properties={})]
        pbf = encode({'name': 'foo', 'features': features},
                     on_invalid_geometry=on_invalid_geometry_make_valid)
        result = decode(pbf)
        features = result['foo']['features']
        self.assertEqual(1, len(features))
        geom = shapely.geometry.shape(features[0]['geometry'])
        self.assertEqual(features[0]['geometry']['type'], 'MultiPolygon')
        self.assertEqual(geom.geom_type, 'MultiPolygon')
        self.assertTrue(geom.is_valid)
        for poly in geom.geoms:
            self.assertTrue(poly.is_valid) 
Example #9
Source File: topology.py    From quantized-mesh-tile with MIT License 6 votes vote down vote up
def _loadGeometry(self, geometrySpec):
        """
        A private method to convert a (E)WKB or (E)WKT to a Shapely geometry.
        """
        if type(geometrySpec) is str and geometrySpec.startswith('POLYGON Z'):
            try:
                geometry = load_wkt(geometrySpec)
            except Exception:
                geometry = None
        else:
            try:
                geometry = load_wkb(geometrySpec)
            except Exception:
                geometry = None

        if geometry is None:
            raise ValueError('Failed to convert WKT or WKB to a Shapely geometry')

        return geometry 
Example #10
Source File: Matchup_test.py    From incubator-sdap-nexus with Apache License 2.0 6 votes vote down vote up
def test_is_pickleable(self):
        edge_point = json.loads("""{
"id": "argo-profiles-5903995(46, 0)",
"time": "2012-10-15T14:24:04Z",
"point": "-33.467 29.728",
"sea_water_temperature": 24.5629997253,
"sea_water_temperature_depth": 2.9796258642,
"wind_speed": null,
"sea_water_salinity": null,
"sea_water_salinity_depth": null,
"platform": 4,
"device": 3,
"fileurl": "ftp://podaac-ftp.jpl.nasa.gov/allData/insitu/L2/spurs1/argo/argo-profiles-5903995.nc"
}""")
        point = DomsPoint.from_edge_point(edge_point)
        self.assertIsNotNone(pickle.dumps(point)) 
Example #11
Source File: test_polygon.py    From mapbox-vector-tile with MIT License 6 votes vote down vote up
def test_polygon_disconnected_inner(self):
        geom = wkt.loads("""POLYGON(
          (0 0, 5 0, 5 5, 0 5, 0 0),
          (1 1, 1 2, 2 2, 1 1),
          (2 1, 2 2, 3 2, 2 1),
          (3 1, 3 2, 4 2, 3 1),
          (1 2, 1 3, 2 3, 1 2),
          (2 2, 2 3, 3 3, 2 2),
          (3 2, 3 3, 4 3, 3 2),
          (1 3, 1 4, 2 4, 1 3),
          (2 3, 2 4, 3 4, 2 3),
          (3 3, 3 4, 4 4, 3 3)
        )""")
        self.assertFalse(geom.is_valid)
        fixed = make_it_valid(geom)
        self.assertTrue(fixed.is_valid)
        self.assertEquals(20.5, fixed.area) 
Example #12
Source File: scihub.py    From OpenSarToolkit with MIT License 6 votes vote down vote up
def create_aoi_str(aoi_wkt):
    '''A helper function to create a scihub API compliant AOI string

    Args:
        aoi (str): is WKT representation of the Area Of Interest

    Returns:
        str: Copernicus' scihub compliant AOI string

    '''

    geom = loads(aoi_wkt)
    if geom.geom_type == 'Point':
        aoi_str = "( footprint:\"Intersects({}, {})\")".format(geom.y, geom.x)

    else:
        # simplify geometry
        aoi_convex = geom.convex_hull

        # create scihub-confrom aoi string
        aoi_str = '( footprint:\"Intersects({})\")'.format(aoi_convex)

    return aoi_str 
Example #13
Source File: vector.py    From OpenSarToolkit with MIT License 6 votes vote down vote up
def gdf_to_json_geometry(gdf):
    """Function to parse features from GeoDataFrame in such a manner 
       that rasterio wants them"""
#    
#    try:
#        gdf.geometry.values[0].type
#        features = [json.loads(gdf.to_json())['features'][0]['geometry']]
#    except AttributeError:
#        ids, feats =[], []
#        for i, feat in enumerate(gdf.geometry.values[0]):
#            ids.append(i)
#            feats.append(feat)
#
#        gdf = gpd.GeoDataFrame({'id': ids,
#                                'geometry': feats}, 
#                                    geometry='geometry', 
#                                    crs = gdf.crs
#                                    )
    geojson = json.loads(gdf.to_json())
    return [feature['geometry'] for feature in geojson['features'] 
            if feature['geometry']] 
Example #14
Source File: test_polygon.py    From mapbox-vector-tile with MIT License 6 votes vote down vote up
def test_polygon_inners_crossing_outer(self):
        geom = wkt.loads("""POLYGON (
          (2325 1015, 2329 1021, 2419 1057, 2461 944, 2369 907, 2325 1015),
          (2329 1012, 2370 909, 2457 944, 2417 1050, 2329 1012),
          (2410 1053, 2410 1052, 2412 1053, 2411 1054, 2410 1053),
          (2378 1040, 2378 1039, 2379 1040, 2379 1041, 2378 1040),
          (2369 1037, 2370 1036, 2371 1036, 2371 1038, 2369 1037),
          (2361 1034, 2362 1033, 2363 1033, 2363 1034, 2361 1034),
          (2353 1031, 2354 1029, 2355 1030, 2354 1031, 2353 1031),
          (2337 1024, 2338 1023, 2339 1023, 2338 1025, 2337 1024)
        )""")
        self.assertFalse(geom.is_valid)
        fixed = make_it_valid(geom)
        self.assertTrue(fixed.is_valid)
        # different versions of GEOS hit this bug in slightly different ways,
        # meaning that some inners get included and some don't, depending on
        # the version. therefore, we need quite a wide range of acceptable
        # answers.
        #
        # the main part of this polygon (outer - largest inner) has area 1551,
        # and the smaller inners sum up to area 11, so we'll take +/-6 from
        # 1545.
        self.assertAlmostEqual(1545, fixed.area, delta=6) 
Example #15
Source File: vector.py    From OpenSarToolkit with MIT License 6 votes vote down vote up
def latlon_to_shp(lon, lat, shapefile):

    shapefile = str(shapefile)

    schema = {'geometry': 'Point',
              'properties': {'id': 'str'}}

    wkt = loads('POINT ({} {})'.format(lon, lat))

    with collection(shapefile, "w",
                    crs=from_epsg(4326),
                    driver="ESRI Shapefile",
                    schema=schema) as output:

        output.write({'geometry': mapping(wkt),
                      'properties': {'id': '1'}}) 
Example #16
Source File: config.py    From mapchete with MIT License 5 votes vote down vote up
def bounds_from_opts(
    wkt_geometry=None, point=None, bounds=None, zoom=None, raw_conf=None
):
    """
    Loads the process pyramid of a raw configuration.

    Parameters
    ----------
    raw_conf : dict
        Raw mapchete configuration as dictionary.

    Returns
    -------
    BufferedTilePyramid
    """
    if wkt_geometry:
        return Bounds(*wkt.loads(wkt_geometry).bounds)
    elif point:
        x, y = point
        zoom_levels = get_zoom_levels(
            process_zoom_levels=raw_conf["zoom_levels"],
            init_zoom_levels=zoom
        )
        tp = raw_conf_process_pyramid(raw_conf)
        return Bounds(*tp.tile_from_xy(x, y, max(zoom_levels)).bounds)
    else:
        return validate_bounds(bounds) if bounds is not None else bounds 
Example #17
Source File: util.py    From geosnap with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def _deserialize_wkb(str):
    return wkb.loads(str, hex=True) 
Example #18
Source File: test_geo.py    From solaris with Apache License 2.0 5 votes vote down vote up
def test_no_overlap(self):
        """Test creation of an overlap object with no intersection."""
        poly_df = pd.read_csv(os.path.join(data_dir, 'sample.csv'))
        polygons = poly_df['PolygonWKT_Pix'].apply(loads).values
        preds = geometries_internal_intersection(polygons)
        # there's no overlap, so result should be an empty GeometryCollection
        assert preds == shapely.geometry.collection.GeometryCollection() 
Example #19
Source File: test_geo.py    From solaris with Apache License 2.0 5 votes vote down vote up
def test_with_overlap(self):
        poly_df = pd.read_csv(os.path.join(data_dir, 'sample.csv'))
        # expand the polygons to generate some overlap
        polygons = poly_df['PolygonWKT_Pix'].apply(
            lambda x: loads(x).buffer(15)).values
        preds = geometries_internal_intersection(polygons)
        with open(os.path.join(data_dir, 'test_overlap_output.txt'), 'r') as f:
            truth = f.read()
            f.close()
        truth = loads(truth)
        # set a threshold for how good overlap with truth has to be in case of
        # rounding errors
        assert truth.intersection(preds).area/truth.area > 0.99 
Example #20
Source File: test_cli.py    From mapchete with MIT License 5 votes vote down vote up
def test_execute_point(mp_tmpdir, example_mapchete, wkt_geom):
    """Using bounds from WKT."""
    g = wkt.loads(wkt_geom)
    run_cli([
        "execute", example_mapchete.path,
        "--point", str(g.centroid.x), str(g.centroid.y)
    ]) 
Example #21
Source File: test_geo.py    From solaris with Apache License 2.0 5 votes vote down vote up
def test_reproject_from_wkt(self):
        input_str = "POLYGON ((736687.5456353347 3722455.06780279, 736686.9301210654 3722464.96326352, 736691.6397869177 3722470.9059681, 736705.5443059544 3722472.614050498, 736706.8992101226 3722462.858909504, 736704.866059878 3722459.457111885, 736713.1443474176 3722452.103498172, 736710.0312805283 3722447.309985571, 736700.3886167214 3722454.263705271, 736698.4577440721 3722451.98534527, 736690.1272768064 3722451.291527834, 736689.4108667439 3722455.113813923, 736687.5456353347 3722455.06780279))"
        result_str = "POLYGON ((-84.4487639 33.6156071, -84.44876790000001 33.6156964, -84.4487156 33.61574889999999, -84.44856540000001 33.6157612, -84.44855339999999 33.61567300000001, -84.44857620000001 33.6156428, -84.448489 33.6155747, -84.4485238 33.6155322, -84.4486258 33.615597, -84.4486472 33.61557689999999, -84.4487371 33.6155725, -84.4487438 33.6156071, -84.4487639 33.6156071))"
        result_geom = loads(result_str)
        reproj_geom = reproject_geometry(input_str, input_crs=32616,
                                         target_crs=4326)
        area_sim = result_geom.intersection(reproj_geom).area/result_geom.area

        assert area_sim > 0.99999 
Example #22
Source File: test_geo.py    From solaris with Apache License 2.0 5 votes vote down vote up
def test_reproject_from_wkt_to_utm(self):
        result_str = "POLYGON ((736687.5456353347 3722455.06780279, 736686.9301210654 3722464.96326352, 736691.6397869177 3722470.9059681, 736705.5443059544 3722472.614050498, 736706.8992101226 3722462.858909504, 736704.866059878 3722459.457111885, 736713.1443474176 3722452.103498172, 736710.0312805283 3722447.309985571, 736700.3886167214 3722454.263705271, 736698.4577440721 3722451.98534527, 736690.1272768064 3722451.291527834, 736689.4108667439 3722455.113813923, 736687.5456353347 3722455.06780279))"
        input_str = "POLYGON ((-84.4487639 33.6156071, -84.44876790000001 33.6156964, -84.4487156 33.61574889999999, -84.44856540000001 33.6157612, -84.44855339999999 33.61567300000001, -84.44857620000001 33.6156428, -84.448489 33.6155747, -84.4485238 33.6155322, -84.4486258 33.615597, -84.4486472 33.61557689999999, -84.4487371 33.6155725, -84.4487438 33.6156071, -84.4487639 33.6156071))"
        result_geom = loads(result_str)
        reproj_geom = reproject_geometry(input_str, input_crs=4326,
                                         target_crs=32616)
        area_sim = result_geom.intersection(reproj_geom).area/result_geom.area

        assert area_sim > 0.99999 
Example #23
Source File: vector.py    From OpenSarToolkit with MIT License 5 votes vote down vote up
def latlon_to_wkt(lat, lon, buffer_degree=None, buffer_meter=None, envelope=False):
    '''A helper function to create a WKT representation of Lat/Lon pair

    This function takes lat and lon vale and returns the WKT Point
    representation by default.

    A buffer can be set in metres, which returns a WKT POLYGON. If envelope
    is set to True, the buffer will be squared by the extent buffer radius.

    Args:
        lat (str): Latitude (deg) of a point
        lon (str): Longitude (deg) of a point
        buffer (float): optional buffer around the point
        envelope (bool): gives a square instead of a circular buffer
                         (only applies if bufferis set)

    Returns:
        wkt (str): WKT string

    '''

    if buffer_degree is None and buffer_meter is None:
        aoi_wkt = 'POINT ({} {})'.format(lon, lat)

    elif buffer_degree:
        aoi_geom = loads('POINT ({} {})'.format(lon, lat)).buffer(buffer_degree)
        if envelope:
            aoi_geom = aoi_geom.envelope

        aoi_wkt = aoi_geom.to_wkt()

    elif buffer_meter:
        aoi_wkt = geodesic_point_buffer(lat, lon, buffer_meter, envelope)

    return aoi_wkt 
Example #24
Source File: create_patches_all.py    From SpaceNet_Off_Nadir_Solutions with Apache License 2.0 5 votes vote down vote up
def patches_and_cocoann(gt, outdir_rgb, outdir_mpan, count=1,create_anns=True):
    gt.index = gt.ImageId
    gt_flt = gt[gt.PolygonWKT_Geo != 'POLYGON EMPTY']
    gv = gt_flt.ImageId.value_counts()
    annotations = []
    images = []
    counter = count
    for u in gt_flt.ImageId.unique():
        try:
            pan_sharpen_dir = ''.join([RAW_DIR, gt.name, '/Pan-Sharpen/'])
            image_file = ''.join([pan_sharpen_dir, 'Pan-Sharpen', '_', u, '.tif'])
            img_rgb = tifffile.imread(image_file)
            if np.argmin(img_rgb.shape) == 0:
                img_rgb = np.transpose(img_rgb, (1, 2, 0))
            img_rgb = img_rgb[:, :, [3, 2, 1, 0]]
            img_mpan = get_pan_sharpend(''.join([RAW_DIR, gt.name, '/']), u)

        except:
            print('load error..', u)
            continue
        if gv[u] > 1:
            poly = gt.loc[u].PolygonWKT_Pix.apply(lambda x: loads(x)).values.tolist()
        else:
            poly = [loads(gt.loc[u].PolygonWKT_Pix)]
        mask = util.mask_for_polygons(poly, im_size=imsize)
        img_patches_rgb, mask_patches, kp = patch_creator.create(img=img_rgb, mask=mask, nonzero=True)
        img_patches_mpan, _, _ = patch_creator.create(img=img_mpan, mask=mask, nonzero=True)

        for i, k in enumerate(kp.keys()):
            file_name_rgb = os.path.join(outdir_rgb, ''.join([u, '_', str(k), '.tif']))
            file_name_mpan = os.path.join(outdir_mpan, ''.join([u, '_', str(k), '.tif']))
            if create_anns:
                anns, images_d, counter = create_coco_anns(file_name_rgb, counter, mask_patches[i].squeeze())
                annotations.extend(anns)
                images.extend(images_d)
            tifffile.imsave(file_name_mpan, img_patches_mpan[i].astype('uint16'))
            tifffile.imsave(file_name_rgb, img_patches_rgb[i].astype('uint16'))
        if DEBUG:
            break
    return annotations, images, counter 
Example #25
Source File: encoder.py    From go2mapillary with GNU General Public License v3.0 5 votes vote down vote up
def _load_geometry(self, geometry_spec):
        if isinstance(geometry_spec, BaseGeometry):
            return geometry_spec

        if isinstance(geometry_spec, dict):
            return SimpleShape(geometry_spec['coordinates'],
                               geometry_spec["type"])

        try:
            return load_wkb(geometry_spec)
        except:
            try:
                return load_wkt(geometry_spec)
            except:
                return None 
Example #26
Source File: idaho.py    From gbdxtools with MIT License 5 votes vote down vote up
def get_images_by_catid(self, catid):
        """ Retrieves the IDAHO image records associated with a given catid.
        Args:
            catid (str): The source catalog ID from the platform catalog.
        Returns:
            results (json): The full catalog-search response for IDAHO images
                            within the catID.
        """

        self.logger.debug('Retrieving IDAHO metadata')

        # get the footprint of the catid's strip
        footprint = self.catalog.get_strip_footprint_wkt(catid)

        # try to convert from multipolygon to polygon:
        try:
            footprint = from_wkt(footprint).geoms[0].wkt
        except:
            pass

        if not footprint:
            self.logger.debug("""Cannot get IDAHO metadata for strip %s,
                                 footprint not found""" % catid)
            return None

        return self.get_images_by_catid_and_aoi(catid=catid,
                                                aoi_wkt=footprint) 
Example #27
Source File: vectors.py    From gbdxtools with MIT License 5 votes vote down vote up
def aggregate_query(self, searchAreaWkt, agg_def, query=None, start_date=None, end_date=None, count=10, index=default_index):
        """Aggregates results of a query into buckets defined by the 'agg_def' parameter.  The aggregations are
        represented by dicts containing a 'name' key and a 'terms' key holding a list of the aggregation buckets.
        Each bucket element is a dict containing a 'term' key containing the term used for this bucket, a 'count' key
        containing the count of items that match this bucket, and an 'aggregations' key containing any child
        aggregations.

        Args:
            searchAreaWkt (str): wkt representation of the geometry
            agg_def (str or AggregationDef): the aggregation definitions
            query (str): a valid Elasticsearch query string to constrain the items going into the aggregation
            start_date (str): either an ISO-8601 date string or a 'now' expression (e.g. "now-6d" or just "now")
            end_date (str): either an ISO-8601 date string or a 'now' expression (e.g. "now-6d" or just "now")
            count (int): the number of buckets to include in the aggregations (the top N will be returned)
            index (str): the index (or alias or wildcard index expression) to run aggregations against, set to None for the entire set of vector indexes

        Returns:
            results (list): A (usually single-element) list of dict objects containing the aggregation results.
        """

        geojson = load_wkt(searchAreaWkt).__geo_interface__
        aggs_str = str(agg_def) # could be string or AggregationDef

        params = {
            "count": count,
            "aggs": aggs_str
        }

        if query:
            params['query'] = query
        if start_date:
            params['start_date'] = start_date
        if end_date:
            params['end_date'] = end_date

        url = self.aggregations_by_index_url % index if index else self.aggregations_url

        r = self.gbdx_connection.post(url, params=params, json=geojson)
        r.raise_for_status()

        return r.json(object_pairs_hook=OrderedDict)['aggregations'] 
Example #28
Source File: vectors.py    From gbdxtools with MIT License 5 votes vote down vote up
def create_from_wkt(self, wkt, item_type, ingest_source, index=None, **attributes):
        '''
        Create a single vector in the vector service

        Args:
            wkt (str): wkt representation of the geometry
            item_type (str): item_type of the vector
            ingest_source (str): source of the vector
            attributes: a set of key-value pairs of attributes
            index (str): optional index to write to, defaults to 'vector-user-provided'

        Returns:
            (str): feature ID
        '''
        # verify the "depth" of the attributes is single layer

        geojson = load_wkt(wkt).__geo_interface__
        vector = {
            'type': "Feature",
            'geometry': geojson,
            'properties': {
                'item_type': item_type,
                'ingest_source': ingest_source,
                'attributes': attributes
            }
        }

        results =  self.create(vector, index=index)
        if len(results['errorMessages']) == 0:
            item = results['successfulItemIds'][0]
            return item.split('/')[-1]
        else:
            raise Exception(results['errorMessages'][0]) 
Example #29
Source File: test_polygon.py    From solaris with Apache License 2.0 5 votes vote down vote up
def test_transform_using_raster(self):
        input_df = os.path.join(data_dir, 'sample.csv')
        input_im = os.path.join(data_dir, 'sample_geotiff.tif')
        output_gdf = georegister_px_df(input_df, im_path=input_im,
                                       geom_col='PolygonWKT_Pix', precision=0)
        truth_df = pd.read_csv(os.path.join(data_dir, 'aff_gdf_result.csv'))
        truth_df['geometry'] = truth_df['geometry'].apply(loads)
        truth_gdf = gpd.GeoDataFrame(
            truth_df,
            crs=rasterio.open(os.path.join(data_dir, 'sample_geotiff.tif')).crs
            )

        assert truth_gdf.equals(output_gdf) 
Example #30
Source File: Matchup_test.py    From incubator-sdap-nexus with Apache License 2.0 5 votes vote down vote up
def test_mur_match(self):
        from shapely.wkt import loads
        from nexustiles.nexustiles import NexusTileService

        polygon = loads("POLYGON((-34.98 29.54, -30.1 29.54, -30.1 31.00, -34.98 31.00, -34.98 29.54))")
        primary_ds = "JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1"
        matchup_ds = "spurs"
        parameter = "sst"
        start_time = 1350259200  # 2012-10-15T00:00:00Z
        end_time = 1350345600  # 2012-10-16T00:00:00Z
        time_tolerance = 86400
        depth_tolerance = 5.0
        radius_tolerance = 1500.0
        platforms = "1,2,3,4,5,6,7,8,9"

        tile_service = NexusTileService()
        tile_ids = [tile.tile_id for tile in
                    tile_service.find_tiles_in_polygon(polygon, primary_ds, start_time, end_time, fetch_data=False,
                                                       fl='id')]
        result = spark_matchup_driver(tile_ids, wkt.dumps(polygon), primary_ds, matchup_ds, parameter, time_tolerance,
                                      depth_tolerance, radius_tolerance, platforms)
        for k, v in result.iteritems():
            print "primary: %s\n\tmatches:\n\t\t%s" % (
                "lon: %s, lat: %s, time: %s, sst: %s" % (k.longitude, k.latitude, k.time, k.sst),
                '\n\t\t'.join(
                    ["lon: %s, lat: %s, time: %s, sst: %s" % (i.longitude, i.latitude, i.time, i.sst) for i in v]))