Python geojson.Polygon() Examples

def way(self, w):
        if not is_polygon(w):
            return

        if "building" not in w.tags:
            return

        if w.tags["building"] in self.building_filter:
            return

        if "location" in w.tags and w.tags["location"] in self.location_filter:
            return

        geometry = geojson.Polygon([[(n.lon, n.lat) for n in w.nodes]])
        shape = shapely.geometry.shape(geometry)

        if shape.is_valid:
            feature = geojson.Feature(geometry=geometry)
            self.storage.add(feature)
        else:
            print("Warning: invalid feature: https://www.openstreetmap.org/way/{}".format(w.id), file=sys.stderr) 

def way(self, w):
        if not is_polygon(w):
            return

        if "amenity" not in w.tags or w.tags["amenity"] != "parking":
            return

        if "parking" in w.tags:
            if w.tags["parking"] in self.parking_filter:
                return

        geometry = geojson.Polygon([[(n.lon, n.lat) for n in w.nodes]])
        shape = shapely.geometry.shape(geometry)

        if shape.is_valid:
            feature = geojson.Feature(geometry=geometry)
            self.storage.add(feature)
        else:
            print("Warning: invalid feature: https://www.openstreetmap.org/way/{}".format(w.id), file=sys.stderr) 

def test_entity_add_city_shape(air_quality_observed):
    air_quality_observed.pop('location')

    air_quality_observed['address']['value'] = {
        "addressCountry": "MX",
        "addressLocality": "Ciudad de México",
    }

    r = geocoding.add_location(air_quality_observed)
    assert r is air_quality_observed

    assert 'location' in r
    assert r['location']['type'] == 'geo:json'

    geo = r['location']['value']
    assert geo['type'] == 'Polygon'
    polygon = geojson.Polygon(geo['coordinates'])
    assert polygon.is_valid 

def crop_aoi(geotiff, aoi, z=0):
    """
    Crop a geographic AOI in a georeferenced image using its RPC functions.

    Args:
        geotiff (string): path or url to the input GeoTIFF image file
        aoi (geojson.Polygon): GeoJSON polygon representing the AOI
        z (float, optional): base altitude with respect to WGS84 ellipsoid (0
            by default)

    Return:
        crop (array): numpy array containing the cropped image
        x, y, w, h (ints): image coordinates of the crop. x, y are the
            coordinates of the top-left corner, while w, h are the dimensions
            of the crop.
    """
    x, y, w, h = bounding_box_of_projected_aoi(rpcm.rpc_from_geotiff(geotiff), aoi, z)
    with rasterio.open(geotiff) as src:
        crop = rasterio_window_crop(src, x, y, w, h)
    return crop, x, y 

def bounding_box_of_projected_aoi(rpc, aoi, z=0, homography=None):
    """
    Return the x, y, w, h pixel bounding box of a projected AOI.

    Args:
        rpc (rpcm.RPCModel): RPC camera model
        aoi (geojson.Polygon): GeoJSON polygon representing the AOI
        z (float): altitude of the AOI with respect to the WGS84 ellipsoid
        homography (2D array, optional): matrix of shape (3, 3) representing an
            homography to be applied to the projected points before computing
            their bounding box.

    Return:
        x, y (ints): pixel coordinates of the top-left corner of the bounding box
        w, h (ints): pixel dimensions of the bounding box
    """
    lons, lats = np.array(aoi['coordinates'][0]).T
    x, y = rpc.projection(lons, lats, z)
    pts = list(zip(x, y))
    if homography is not None:
        pts = points_apply_homography(homography, pts)
    return np.round(bounding_box2D(pts)).astype(int) 

def geojson_lonlat_to_utm(aoi):
    """
    """
    # compute the utm zone number of the first polygon vertex
    lon, lat = aoi['coordinates'][0][0]
    utm_zone = utm.from_latlon(lat, lon)[2]

    # convert all polygon vertices coordinates from (lon, lat) to utm
    c = []
    for lon, lat in aoi['coordinates'][0]:
        c.append(utm.from_latlon(lat, lon, force_zone_number=utm_zone)[:2])

    return geojson.Polygon([c]) 

def valid_geojson(filepath):
    """
    Check if a file contains valid geojson.
    """
    with open(filepath, 'r') as f:
        geo = geojson.load(f)
    if type(geo) == geojson.geometry.Polygon:
        return geo
    if type(geo) == geojson.feature.Feature:
        p = geo['geometry']
        if type(p) == geojson.geometry.Polygon:
            return p
    if type(geo) == geojson.feature.FeatureCollection:
        p = geo['features'][0]['geometry']
        if type(p) == geojson.geometry.Polygon:
            return p
    raise argparse.ArgumentTypeError('Invalid geojson: only polygons are supported') 

def geometry(self):
        try:
            if self.type() == 'node':
                if not self.lon() or not self.lat():
                    self._raiseException('Cannot build geometry: geometry information not included.')
                return geojson.Point((self.lon(), self.lat()))
            elif self.type() == 'way':
                if not self.__getElement('geometry'):
                    self._raiseException('Cannot build geometry: geometry information not included.')
                cs = self.__geometry_csToList(self.__getElement('geometry'))
                if self.__geometry_equal(cs[0], cs[-1]):
                    return geojson.Polygon([cs])
                else:
                    return geojson.LineString(cs)
            elif self.type() == 'relation':
                members = copy.deepcopy(self.__members())
                membersOuter = self.__geometry_extract(members, 'outer')
                if len(membersOuter) == 0:
                    self._raiseException('Cannot build geometry: no outer rings found.')
                membersInner = self.__geometry_extract(members, 'inner')
                ringsOuter = self.__geometry_buildRings(membersOuter)
                ringsInner = self.__geometry_buildRings(membersInner)
                ringsOuter = self.__geometry_orientRings(ringsOuter, positive=True)
                ringsInner = self.__geometry_orientRings(ringsInner, positive=False)
                polygons = self.__geometry_buildPolygons(ringsOuter, ringsInner)
                if len(polygons) > 1:
                    return geojson.MultiPolygon(polygons)
                else:
                    return geojson.Polygon(polygons[0])
            else:
                self._raiseException('Cannot build geometry: type of element unknown.')
        except Exception as e:
            _extendAndRaiseException(e, ' ({}/{})'.format(self.type(), self.id())) 

def coords_series_to_geometry(coords, geomtype='linestring', dtype='geojson'):
    """
    Convert a series of coords (list of list(s)) to a series of geometry objects.
    :param coords: series of lists of xy coordinates
    :param geomtype: geometry type of target 
    :param dtype: format of geometry objects to be created ('geojson', 'shapely')
    :return: series of geometry objects
    >>> import swmmio
    >>> model = swmmio.Model(MODEL_FULL_FEATURES_XY)
    >>> nodes = model.nodes()
    >>> geoms = coords_series_to_geometry(nodes['coords'], geomtype='point')
    >>> geoms.iloc[0]
    {"coordinates": [2748073.3060000003, 1117746.087], "type": "Point"}
    """

    # detect whether LineString or Point should be used
    geomtype = geomtype.lower()
    if geomtype == 'linestring':
        geoms = [LineString(latlngs) for latlngs in coords]
    elif geomtype == 'point':
        geoms = [Point(latlngs[0]) for latlngs in coords]
    elif geomtype == 'polygon':
        geoms = [Polygon([latlngs]) for latlngs in coords]

    if dtype.lower() == 'shape':
        # convert to shapely objects
        try:
            from shapely.geometry import shape
        except ImportError:
            raise ImportError('shapely module needed. Install it via GeoPandas with conda: ',
                              'conda install geopandas')
        geoms = [shape(g) for g in geoms]

    return pd.Series(index=coords.index, name='geometry', data=geoms) 

def utm_bbx(aoi, epsg=None, r=None, offset=(0, 0)):
    """
    Compute UTM bounding box of a longitude, latitude AOI.

    Args:
        aoi (geojson.Polygon): area of interest, defined as a (lon, lat) polygon
        epsg (int): EPSG code of the desired UTM zone
        r (int): if not None, round bounding box vertices to vertices of an
            r-periodic grid
        offset (tuple): origin of the r-periodic grid

    Returns:
        ulx, uly, lrx, lry (floats): bounding box upper left (ul) and lower
            right (lr) x, y coordinates
        epsg (int): EPSG code of the UTM zone
    """
    if epsg is None:  # compute the EPSG code of the AOI centroid
        lon, lat = np.mean(aoi['coordinates'][0][:-1], axis=0)
        epsg = compute_epsg(lon, lat)

    # convert all polygon vertices coordinates from (lon, lat) to utm
    lons, lats = np.asarray(aoi['coordinates'][0]).T
    xs, ys = pyproj_transform(lons, lats, 4326, epsg)
    c = list(zip(xs, ys))

    # utm bounding box
    x, y, w, h = bounding_box2D(c)  # minx, miny, width, height
    ulx, uly, lrx, lry = x, y + h, x + w, y

    if r is not None:  # round to multiples of the given resolution
        ox, oy = offset
        ulx = ox + r * np.round((ulx - ox) / r)
        uly = oy + r * np.round((uly - oy) / r)
        lrx = ox + r * np.round((lrx - ox) / r)
        lry = oy + r * np.round((lry - oy) / r)

    return ulx, uly, lrx, lry, epsg 

def geojson_geometry_object(lat, lon, w, h):
    """
    """
    return geojson.Polygon([lonlat_rectangle_centered_at(lon, lat, w, h)]) 

def download(imgs, bands, aoi, mirror, out_dir, parallel_downloads):
    """
    Download a timeseries of crops with GDAL VSI feature.

    Args:
        imgs (list): list of images
        bands (list): list of bands
        aoi (geojson.Polygon): area of interest
        mirror (str): either 'aws' or 'gcloud'
        out_dir (str): path where to store the downloaded crops
        parallel_downloads (int): number of parallel downloads
    """
    coords = utils.utm_bbx(aoi)
    crops_args = []
    for img in imgs:
        for b in set(bands + ['BQA']):
            fname = os.path.join(out_dir, '{}_band_{}.tif'.format(img.filename, b))
            crops_args.append((fname, img.urls[mirror][b]))

    out_dir = os.path.abspath(os.path.expanduser(out_dir))
    os.makedirs(out_dir, exist_ok=True)
    print('Downloading {} crops ({} images with {} bands)...'.format(len(crops_args),
                                                                     len(imgs),
                                                                     len(bands) +1),
          end=' ')
    parallel.run_calls(utils.rasterio_geo_crop, crops_args,
                       extra_args=(*coords,), pool_type='threads',
                       nb_workers=parallel_downloads) 

def search(aoi, start_date=None, end_date=None, satellite='L8',
           sensor='OLITIRS', api='devseed'):
    """
    Search Landsat images covering an AOI and timespan using a given API.

    Args:
        aoi (geojson.Polygon): area of interest
        start_date (datetime.datetime): start of the search time range
        end_date (datetime.datetime): end of the search time range
        api (str, optional): either gcloud (default), scihub, planet or devseed
        satellite (str, optional): either L1...L8
        sensor (str, optional): MSS, TM, ETM, OLITIRS
        see https://landsat.usgs.gov/what-are-band-designations-landsat-satellites

    Returns:
        list of image objects
    """
    # list available images
    if api == 'gcloud':
        from tsd import search_gcloud
        images = search_gcloud.search(aoi, start_date, end_date, satellite=satellite, sensor=sensor)
    elif api == 'devseed':
        from tsd import search_devseed
        images = search_devseed.search(aoi, start_date, end_date, satellite='Landsat-8')
    else:
        raise ValueError('The api "{}" is not available for {}.'.format(api, __file__))

    images = [l8_metadata_parser.LandsatImage(img, api) for img in images]

    # sort images by acquisition day, then by mgrs id
    images.sort(key=(lambda k: (k.date.date(), k.row, k.path)))

    print('Found {} images'.format(len(images)))
    return images 

def write_to(data, property_names, output_file):
    '''
    Write list of tuples to geojson.
       First entry of each tuple should be geometry in hex coordinates
       and the rest properties.

       Args:
           data: List of tuples.
           property_names: List of strings. Should be same length as the
                           number of properties.
           output_file (str): Output file name.

    '''

    geojson_features = []
    for entry in data:
        coords_in_hex, properties = entry[0], entry[1:]
        geometry = loads(coords_in_hex, hex=True)
        property_dict = dict(zip(property_names, properties))
        if geometry.geom_type == 'Polygon':
            coords = [list(geometry.exterior.coords)]   # brackets required
            geojson_feature = geojson.Feature(geometry=geojson.Polygon(coords),
                                              properties=property_dict)
        elif geometry.geom_type == 'Point':
            coords = list(geometry.coords)[0]
            geojson_feature = geojson.Feature(geometry=geojson.Point(coords),
                                              properties=property_dict)
        geojson_features.append(geojson_feature)

    feature_collection = geojson.FeatureCollection(geojson_features)

    with open(output_file, 'wb') as f:
        geojson.dump(feature_collection, f) 

def test_compile_in_op(self):
        """Test private function _compile_in_op"""
        self.assertEqual(_compile_in_op('a', ['a', 'b']), 'p.get("a") in [\'a\', \'b\']')
        self.assertEqual(_compile_in_op('$type', ['Point', 'Polygon']), 'f.get("geometry").get("type") in [\'Point\', \'Polygon\']') 

def test_compile_comparison_op(self):
        """Test private function _compile_comparison_op"""
        self.assertEqual(_compile_comparison_op('a', 5, '=='), 'p.get("a")==5')
        self.assertEqual(_compile_comparison_op('$type', 'Polygon', '=='), 'f.get("geometry").get("type")=="Polygon"') 

def test_geometry_in(self):
        """Test $type specific filters for inclusion"""
        ff = create_filter(['in', '$type', 'Point', 'Polygon'])
        # print(_compile(['in', '$type', 'Point', 'Polygon']))
        passing = Feature(geometry=polygon_geometry)
        failing = Feature(geometry=line_geometry)
        self.assertTrue(ff(passing))
        self.assertFalse(ff(failing)) 

def test_geometry_comparison(self):
        """Test $type specific filters for comparison"""
        ff = create_filter(['==', '$type', 'Polygon'])
        # print(_compile(['==', '$type', 'Polygon']))
        passing = Feature(geometry=polygon_geometry)
        failing = Feature(geometry=line_geometry)
        self.assertTrue(ff(passing))
        self.assertFalse(ff(failing)) 

def load_polygon_netcdf(path):
        """
        Loads polygons from netcdf
        :param path: absolute path and filename
        :return: geojson object for plotting in bokeh
        """

        return getRDT(path, 0, 'Polygon') 

def box_to_json_rep(box: SlfBox) -> Polygon:
    ps = list_point_tuples(box.to_polygon())
    return Polygon(ps) 

def polygon_to_json_rep(polygon: SlfPolygon) -> Polygon:
    ps = list_point_tuples(polygon)
    return Polygon(ps) 

def contourf_to_geojson_overlap(contourf, geojson_filepath=None, min_angle_deg=None,
                                ndigits=5, unit='', stroke_width=1, fill_opacity=.9,
                                geojson_properties=None, strdump=False, serialize=True):
    """Transform matplotlib.contourf to geojson with overlapping filled contours."""
    polygon_features = []
    contourf_idx = 0
    contourf_levels = get_contourf_levels(contourf.levels, contourf.extend)
    for collection in contourf.collections:
        color = collection.get_facecolor()
        for path in collection.get_paths():
            for coord in path.to_polygons():
                if min_angle_deg:
                    coord = keep_high_angle(coord, min_angle_deg)
                coord = np.around(coord, ndigits) if ndigits else coord
                polygon = Polygon(coordinates=[coord.tolist()])
                fcolor = rgb2hex(color[0])
                properties = set_contourf_properties(stroke_width, fcolor, fill_opacity, contourf_levels[contourf_idx], unit)
                if geojson_properties:
                    properties.update(geojson_properties)
                feature = Feature(geometry=polygon, properties=properties)
                polygon_features.append(feature)
        contourf_idx += 1
    feature_collection = FeatureCollection(polygon_features)
    return _render_feature_collection(feature_collection, geojson_filepath, strdump, serialize) 

def get_time_series(aoi, start_date=None, end_date=None, bands=['B8'],
                    satellite='Landsat', sensor=None,
                    out_dir='', api='devseed', mirror='gcloud',
                    cloud_masks=False, check_empty=False,
                    parallel_downloads=multiprocessing.cpu_count()):
    """
    Main function: crop and download a time series of Sentinel-2 images.

    Args:
        aoi (geojson.Polygon): area of interest
        start_date (datetime.datetime, optional): start of the time range
        end_date (datetime.datetime, optional): end of the time range
        bands (list, optional): list of bands
        satellite (str, optional): either L1...L8
        sensor (str, optional): MSS, TM, ETM, OLI
        see https://landsat.usgs.gov/what-are-band-designations-landsat-satellites
        out_dir (str, optional): path where to store the downloaded crops
        api (str, optional): either devseed (default), scihub, planet or gcloud
        mirror (str, optional): either 'aws' or 'gcloud'
        cloud_masks (bool, optional): if True, cloud masks are downloaded and
            cloudy images are discarded
        check_empty (bool, optional): if True, QA masks are downloaded and
            empty images are discarded
        parallel_downloads (int): number of parallel gml files downloads
    """
    # check access to the selected search api and download mirror
    check_args(api, mirror)

    # default date range
    if end_date is None:
        end_date = datetime.date.today()
    if start_date is None:
        start_date = end_date - datetime.timedelta(91)  # 3 months

    # list available images
    images = search(aoi, start_date, end_date, satellite, sensor, api=api)

    # download crops
    download(images, bands, aoi, mirror, out_dir, parallel_downloads)

    # discard images that failed to download
    images = [i for i in images if bands_files_are_valid(i, bands, api, out_dir)]

    # embed all metadata as GeoTIFF tags in the image files
    for img in images:
        img['downloaded_by'] = 'TSD on {}'.format(datetime.datetime.now().isoformat())
        for b in bands:
            filepath = os.path.join(out_dir, '{}_band_{}.tif'.format(img.filename, b))
            utils.set_geotif_metadata_items(filepath, img)

    if cloud_masks:  # discard images that are totally covered by clouds
        read_cloud_masks(images, bands, parallel_downloads,
                         out_dir=out_dir)

    if check_empty:  # discard images that are totally empty
        read_empty_images(images, bands, parallel_downloads,
                         out_dir=out_dir) 

def process_solr_tiles(self, solr_tiles):
        """ processes the solr_json 
            discovery geo tiles,
            aggregating to a certain
            depth
        """
        # first aggregate counts for tile that belong togther
        aggregate_tiles = self.aggregate_spatial_tiles(solr_tiles)
        # now generate GeoJSON for each tile region
        # print('Total tiles: ' + str(t) + ' reduced to ' + str(len(aggregate_tiles)))
        i = 0
        for tile_key, aggregate_count in aggregate_tiles.items():
            i += 1
            add_region = True
            fl = FilterLinks()
            fl.base_request_json = self.filter_request_dict_json
            fl.spatial_context = self.spatial_context
            new_rparams = fl.add_to_request('disc-geotile',
                                            tile_key)
            record = LastUpdatedOrderedDict()
            record['id'] = fl.make_request_url(new_rparams)
            record['json'] = fl.make_request_url(new_rparams, '.json')
            record['count'] = aggregate_count
            record['type'] = 'Feature'
            record['category'] = 'oc-api:geo-facet'
            if self.min_date is not False \
               and self.max_date is not False:
                when = LastUpdatedOrderedDict()
                when['id'] = '#event-' + tile_key
                when['type'] = 'oc-gen:formation-use-life'
                # convert numeric to GeoJSON-LD ISO 8601
                when['start'] = ISOyears().make_iso_from_float(self.min_date)
                when['stop'] = ISOyears().make_iso_from_float(self.max_date)
                record['when'] = when
            gm = GlobalMercator()
            geo_coords = gm.quadtree_to_geojson_poly_coords(tile_key)
            geometry = LastUpdatedOrderedDict()
            geometry['id'] = '#geo-disc-tile-geom-' + tile_key
            geometry['type'] = 'Polygon'
            geometry['coordinates'] = geo_coords
            record['geometry'] = geometry
            properties = LastUpdatedOrderedDict()
            properties['id'] = '#geo-disc-tile-' + tile_key
            properties['href'] = record['id']
            properties['label'] = 'Discovery region (' + str(i) + ')'
            properties['feature-type'] = 'discovery region (facet)'
            properties['count'] = aggregate_count
            properties['early bce/ce'] = self.min_date
            properties['late bce/ce'] = self.max_date
            record['properties'] = properties
            if len(tile_key) >= 6:
                if tile_key[:6] == '211111':
                    # no bad coordinates (off 0, 0 coast of Africa)
                    add_region = False  # don't display items without coordinates
            if add_region:
                self.geojson_regions.append(record) 

def export_geojson(self, filename):
        import geojson
        self._log.debug('Exporting GeoJSON Quadtree to %s', filename)
        features = []

        for lf in self.leaves:
            llN, llE, urN, urE = (lf.llN, lf.llE, lf.urN, lf.urE)

            if self.frame.isDegree():
                llN += self.frame.llLat
                llE += self.frame.llLon
                urN += self.frame.llLat
                urE += self.frame.llLon

            coords = num.array([
                (llN, llE),
                (llN, urE),
                (urN, urE),
                (urN, llE),
                (llN, llE)])

            if self.frame.isMeter():
                coords = od.ne_to_latlon(
                    self.frame.llLat, self.frame.llLon, *coords.T)
                coords = num.array(coords).T

            coords = coords[:, [1, 0]].tolist()

            feature = geojson.Feature(
                geometry=geojson.Polygon(coordinates=[coords]),
                id=lf.id,
                properties={
                    'mean': float(lf.mean),
                    'median': float(lf.median),
                    'std': float(lf.std),
                    'var': float(lf.var),

                    'phi': float(lf.phi),
                    'theta': float(lf.theta),
                    'unitE': float(lf.unitE),
                    'unitN': float(lf.unitN),
                    'unitU': float(lf.unitU),
                })
            features.append(feature)

        collection = geojson.FeatureCollection(
            features)
        with open(filename, 'w') as f:
            geojson.dump(collection, f)