Python shapely.wkt() Examples

def wkt_to_shp(wkt_list, shp_file):
    '''Take output of build_graph_wkt() and render the list of linestrings
    into a shapefile
    # https://gis.stackexchange.com/questions/52705/how-to-write-shapely-geometries-to-shapefiles
    '''

    # Define a linestring feature geometry with one attribute
    schema = {
        'geometry': 'LineString',
        'properties': {'id': 'int'},
    }

    # Write a new shapefile
    with fiona.open(shp_file, 'w', 'ESRI Shapefile', schema) as c:
        for i, line in enumerate(wkt_list):
            shape = shapely.wkt.loads(line)
            c.write({
                    'geometry': mapping(shape),
                    'properties': {'id': i},
                    })

    return


############################################################################### 

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 

def snwe2file(snwe):
    """Use Shapely to convert to GeoJSON & WKT.

    Save local text files in variety of formats to record bounds: snwe.json,
    snwe.wkt, snwe.txt.

    Parameters
    ----------
    snwe : list
        bounding coordinates [south, north, west, east].

    """
    S, N, W, E = snwe
    roi = box(W, S, E, N)
    with open("snwe.json", "w") as j:
        json.dump(mapping(roi), j)
    with open("snwe.wkt", "w") as w:
        w.write(roi.wkt)
    with open("snwe.txt", "w") as t:
        snweList = "[{0:.3f}, {1:.3f}, {2:.3f}, {3:.3f}]".format(S, N, W, E)
        t.write(snweList) 

def getPointCoords(row, geom, coord_type):
    """Calculates coordinates ('x' or 'y') of a Point geometry"""
    geom_val=row[geom]
    if isinstance(geom_val, str):
        P = shapely.wkt.loads(geom_val)
    else:P=geom_val    
    # print(type(P))  #P.is_empty
    # print(P.x)
    # print("+"*50)
    # if isinstance(P,float):
        # print(P)
        # P=Point()
    # print(P)
    if coord_type == 'x':
        if isinstance(P,float): return float('nan')
        else: return P.x
    elif coord_type == 'y':
        if isinstance(P,float): return float('nan')
        else: return P.y 

def to_geodataframe(products):
        """Return the products from a query response as a GeoPandas GeoDataFrame
        with the values in their appropriate Python types.
        """
        try:
            import geopandas as gpd
            import shapely.wkt
        except ImportError:
            raise ImportError(
                "to_geodataframe requires the optional dependencies GeoPandas and Shapely."
            )

        crs = "EPSG:4326"  # WGS84
        if len(products) == 0:
            return gpd.GeoDataFrame(crs=crs, geometry=[])

        df = SentinelAPI.to_dataframe(products)
        geometry = [shapely.wkt.loads(fp) for fp in df["footprint"]]
        # remove useless columns
        df.drop(["footprint", "gmlfootprint"], axis=1, inplace=True)
        return gpd.GeoDataFrame(df, crs=crs, geometry=geometry) 

def to_geojson(products):
        """Return the products from a query response as a GeoJSON with the values in their
        appropriate Python types.
        """
        feature_list = []
        for i, (product_id, props) in enumerate(products.items()):
            props = props.copy()
            props["id"] = product_id
            poly = geomet.wkt.loads(props["footprint"])
            del props["footprint"]
            del props["gmlfootprint"]
            # Fix "'datetime' is not JSON serializable"
            for k, v in props.items():
                if isinstance(v, (date, datetime)):
                    props[k] = v.strftime("%Y-%m-%dT%H:%M:%S.%fZ")
            feature_list.append(geojson.Feature(geometry=poly, id=i, properties=props))
        return geojson.FeatureCollection(feature_list) 

def _parse_geometry(geometry):
        """ Parses given geometry into shapely object

        :param geometry:
        :return: Shapely polygon or multipolygon
        :rtype: shapely.geometry.Polygon or shapely.geometry.MultiPolygon
        :raises TypeError
        """
        if isinstance(geometry, str):
            geometry = shapely.wkt.loads(geometry)
        elif isinstance(geometry, dict):
            geometry = shapely.geometry.shape(geometry)
        elif not isinstance(geometry, shapely.geometry.base.BaseGeometry):
            raise TypeError('Unsupported geometry representation')

        if not isinstance(geometry, (shapely.geometry.Polygon, shapely.geometry.MultiPolygon)):
            raise ValueError('Supported geometry types are polygon and multipolygon, got {}'.format(type(geometry)))

        return geometry 

def convert_pix_lstring_to_geo(wkt_lstring, im_file):
    '''Convert linestring in pixel coords to geo coords'''
    shape = wkt_lstring  # shapely.wkt.loads(lstring)
    x_pixs, y_pixs = shape.coords.xy
    coords_latlon = []
    coords_utm = []
    for (x, y) in zip(x_pixs, y_pixs):
        lon, lat = apls_utils.pixelToGeoCoord(x, y, im_file)
        [utm_east, utm_north, utm_zone, utm_letter] = utm.from_latlon(lat, lon)
        coords_utm.append([utm_east, utm_north])
        coords_latlon.append([lon, lat])

    lstring_latlon = LineString([Point(z) for z in coords_latlon])
    lstring_utm = LineString([Point(z) for z in coords_utm])

    return lstring_latlon, lstring_utm, utm_zone, utm_letter


############################################################################### 

def __repr__(self):
        """ Method for class representation
        """
        return '{}({}, crs={})'.format(self.__class__.__name__, self.wkt, self.crs) 

def wkt(self):
        """ Transforms geometry object into `Well-known text` format

        :return: string in WKT format
        :rtype: str
        """
        return self.geometry.wkt 

def get_edge_geo_coords(G, im_file, remove_pix_geom=True,
                        verbose=False):

    ne = len(list(G.edges()))
    for i, (u, v, attr_dict) in enumerate(G.edges(data=True)):
        if verbose:
            print("edge:", u, v)
        if (i % 1000) == 0:
            print("edge", i, "/", ne)
        geom_pix = attr_dict['geometry_pix']
        lstring_latlon, lstring_utm, utm_zone, utm_letter = convert_pix_lstring_to_geo(
            geom_pix, im_file)
        attr_dict['geometry_latlon_wkt'] = lstring_latlon.wkt
        attr_dict['geometry_utm_wkt'] = lstring_utm.wkt
        attr_dict['length_latlon'] = lstring_latlon.length
        attr_dict['length_utm'] = lstring_utm.length
        attr_dict['length'] = lstring_utm.length
        attr_dict['utm_zone'] = utm_zone
        attr_dict['utm_letter'] = utm_letter
        if verbose:
            print("  attr_dict:", attr_dict)

        # geometry screws up osmnx.simplify function
        if remove_pix_geom:
            #attr_dict['geometry_wkt'] = lstring_latlon.wkt
            attr_dict['geometry_pix'] = geom_pix.wkt

    return G


############################################################################### 

def placename_to_wkt(placename):
    """Geocodes the placename to rectangular bounding extents using Nominatim API and
       returns the corresponding 'ENVELOPE' form Well-Known-Text.

    Parameters
    ----------
    placename : str
        the query to geocode

    Returns
    -------
    full name and coordinates of the queried placename
        list in the form [wkt string in form 'ENVELOPE(minX, maxX, maxY, minY)', string placeinfo]
    """

    rqst = requests.post(
        "https://nominatim.openstreetmap.org/search", params={"q": placename, "format": "geojson"}
    )
    # Check that the response from Openstreetmapserver has status code 2xx
    # and that the response is valid JSON.
    rqst.raise_for_status()
    jsonlist = rqst.json()
    if len(jsonlist["features"]) == 0:
        raise ValueError('Unable to find a matching location for "{}"'.format(placename))
    # Get the First result's bounding box and description.
    feature = jsonlist["features"][0]
    minX, minY, maxX, maxY = feature["bbox"]
    footprint = "ENVELOPE({}, {}, {}, {})".format(minX, maxX, maxY, minY)
    placeinfo = feature["properties"]["display_name"]
    return footprint, placeinfo 

def _internal_test(mask_dir, out_file):

    # Postprocessing phase

    fn_out = out_file
    with open(fn_out, 'w') as f:
        f.write("ImageId,BuildingId,PolygonWKT_Pix,Confidence\n")
        test_image_list = os.listdir(os.path.join(mask_dir))
        for idx, image_id in tqdm(enumerate(test_image_list),
                                       total=len(test_image_list)):
            img1 = cv2.imread(os.path.join(mask_dir, image_id), cv2.IMREAD_UNCHANGED)
            labels = img1.astype(np.uint16)
            df_poly = mask_to_poly(labels, min_polygon_area_th=MIN_AREA)
            if len(df_poly) > 0:
                for i, row in df_poly.iterrows():
                    line = "{},{},\"{}\",{:.6f}\n".format(
                        image_id.lstrip("Pan-Sharpen_").rstrip(".tif"),
                        row.bid,
                        row.wkt,
                        row.area_ratio)
                    line = _remove_interiors(line)
                    f.write(line)
            else:
                f.write("{},{},{},0\n".format(
                    image_id,
                    -1,
                    "POLYGON EMPTY")) 

def mask_to_poly(mask, min_polygon_area_th=MIN_AREA):
    shapes = rasterio.features.shapes(mask.astype(np.int16), mask > 0)
    poly_list = []
    mp = shapely.ops.cascaded_union(
        shapely.geometry.MultiPolygon([
            shapely.geometry.shape(shape)
            for shape, value in shapes
        ]))

    if isinstance(mp, shapely.geometry.Polygon):
        df = pd.DataFrame({
            'area_size': [mp.area],
            'poly': [mp],
        })
    else:
        df = pd.DataFrame({
            'area_size': [p.area for p in mp],
            'poly': [p for p in mp],
        })

    df = df[df.area_size > min_polygon_area_th].sort_values(
        by='area_size', ascending=False)
    df.loc[:, 'wkt'] = df.poly.apply(lambda x: shapely.wkt.dumps(
        x, rounding_precision=0))
    df.loc[:, 'bid'] = list(range(1, len(df) + 1))
    df.loc[:, 'area_ratio'] = df.area_size / df.area_size.max()
    return df 

def format(cls, value: shapely.geometry.base.BaseGeometry) -> str:
        return value.wkt if value else '' 

def format(cls, value: Optional[shapely.geometry.Polygon]) -> str:
        return value.wkt if value else '' 

def load_asf_json(jsonfile):
    """Convert JSON metadata from ASF query to dataframe.

    JSON metadata returned from ASF DAAC API is loaded into a geopandas
    GeoDataFrame, with timestamps converted to datatime objects.

    Parameters
    ----------
    jsonfile : str
        Path to the json file from an ASF API query.

    Returns
    -------
    gf :  GeoDataFrame
        A geopandas GeoDataFrame

    """
    with open(jsonfile) as f:
        meta = json.load(f)[0]  # list of scene dictionaries

    df = pd.DataFrame(meta)
    polygons = df.stringFootprint.apply(shapely.wkt.loads)
    gf = gpd.GeoDataFrame(df, crs="EPSG:4326", geometry=polygons)

    gf["timeStamp"] = pd.to_datetime(gf.sceneDate, format="%Y-%m-%d %H:%M:%S")
    gf["sceneDateString"] = gf.timeStamp.apply(lambda x: x.strftime("%Y-%m-%d"))
    gf["dateStamp"] = pd.to_datetime(gf.sceneDateString)
    gf["utc"] = gf.timeStamp.apply(lambda x: x.strftime("%H:%M:%S"))
    gf["orbitCode"] = gf.relativeOrbit.astype("category").cat.codes

    return gf 

def geojson_to_wkt(geojson_obj, feature_number=0, decimals=4):
    """Convert a GeoJSON object to Well-Known Text. Intended for use with OpenSearch queries.

    In case of FeatureCollection, only one of the features is used (the first by default).
    3D points are converted to 2D.

    Parameters
    ----------
    geojson_obj : dict
        a GeoJSON object
    feature_number : int, optional
        Feature to extract polygon from (in case of MultiPolygon
        FeatureCollection), defaults to first Feature
    decimals : int, optional
        Number of decimal figures after point to round coordinate to. Defaults to 4 (about 10
        meters).

    Returns
    -------
    polygon coordinates
        string of comma separated coordinate tuples (lon, lat) to be used by SentinelAPI
    """
    if "coordinates" in geojson_obj:
        geometry = geojson_obj
    elif "geometry" in geojson_obj:
        geometry = geojson_obj["geometry"]
    else:
        geometry = geojson_obj["features"][feature_number]["geometry"]

    def ensure_2d(geometry):
        if isinstance(geometry[0], (list, tuple)):
            return list(map(ensure_2d, geometry))
        else:
            return geometry[:2]

    def check_bounds(geometry):
        if isinstance(geometry[0], (list, tuple)):
            return list(map(check_bounds, geometry))
        else:
            if geometry[0] > 180 or geometry[0] < -180:
                raise ValueError("Longitude is out of bounds, check your JSON format or data")
            if geometry[1] > 90 or geometry[1] < -90:
                raise ValueError("Latitude is out of bounds, check your JSON format or data")

    # Discard z-coordinate, if it exists
    geometry["coordinates"] = ensure_2d(geometry["coordinates"])
    check_bounds(geometry["coordinates"])

    wkt = geomet.wkt.dumps(geometry, decimals=decimals)
    # Strip unnecessary spaces
    wkt = re.sub(r"(?<!\d) ", "", wkt)
    return wkt 

def to_geometry(cls, value: Any, geom_type: type = shapely.geometry.base.BaseGeometry) \
            -> Optional[shapely.geometry.base.BaseGeometry]:
        if value is None:
            return None

        if isinstance(value, geom_type):
            # Compatible!
            # noinspection PyTypeChecker
            return value

        if isinstance(value, shapely.geometry.base.BaseGeometry):
            # Incompatible!
            raise ValidationError(cls.errmsg("Passed geometry type is incompatible."))

        # Try converting to compatible type
        if isinstance(value, str):
            value = value.strip()
            if value == '':
                # Empty strings are same as None
                return None

            try:
                try:
                    coordinates = [float(v) for v in value.split(',', maxsplit=3)]
                    if len(coordinates) == 2 and cls._is_compatible_type(geom_type, shapely.geometry.Point):
                        # Point "x, y"?
                        x, y = coordinates
                        return shapely.geometry.Point(x, y)
                    elif len(coordinates) == 4 and cls._is_compatible_type(geom_type, shapely.geometry.Polygon):
                        # Polygon box "x1, y1, x2, y2"?
                        x1, y1, x2, y2 = coordinates
                        return shapely.geometry.box(x1, y1, x2, y2)
                except (ValueError, TypeError):
                    # Geometry WKT?
                    return shapely.wkt.loads(value)
            except (ShapelyError, ValueError, TypeError) as e:
                raise ValidationError(cls.errmsg('Invalid geometry WKT format.')) from e

            raise ValidationError(cls.errmsg("Unrecognized text format."))

        if cls._is_compatible_type(geom_type, shapely.geometry.Point):
            try:
                # Coordinates list that forms a polygon?
                return shapely.geometry.Point(value)
            except (ShapelyError, ValueError, TypeError):
                pass

        if cls._is_compatible_type(geom_type, shapely.geometry.Polygon):
            try:
                # Coordinates list that forms a polygon?
                return shapely.geometry.Polygon(value)
            except (ShapelyError, ValueError, TypeError):
                # Polygon box x1, y1, x2, y2?
                try:
                    x1, y1, x2, y2 = value
                    return shapely.geometry.box(x1, y1, x2, y2)
                except (ShapelyError, ValueError, TypeError) as e:
                    raise ValidationError(cls.errmsg(str(e))) from e

        raise ValidationError(cls.errmsg()) 

def covid_19_csv2gpd(dataFpDic):
    Covid_19=pd.read_csv(dataFpDic["Covid-19_CasesByZipCode"])
    covid19_dfCopy=Covid_19.copy(deep=True)
    # print("-"*50)
    # print(Covid_19.head)
    # print(Covid_19.columns)
    # Covid_19_MultiIdx_weekZip=Covid_19
    Covid_19['zipBak']=Covid_19['ZIP Code']
    covid19_df_byZip=covid19_dfCopy.set_index(['ZIP Code','Week Number']).sort_index()
    
    Covid_19_MultiIdx_weekZip=Covid_19.set_index(['Week Number','ZIP Code']).sort_index()
    Covid_19_MultiIdx_weekZip=Covid_19_MultiIdx_weekZip.rename(columns={
    'Week Start':'WeekStart', 
    'Week End':'WeekEnd', 
    'Cases - Weekly':'CasesWeekly', 
    'Cases - Cumulative':'CasesCumulative',
    'Case Rate - Weekly':'CaseRateWeekly', 
    'Case Rate - Cumulative':'CaseRateCumulative', 
    'Tests - Weekly':'TestsWeekly',
    'Tests - Cumulative':'TestsCumulative', 
    'Test Rate - Weekly':'TestRateWeekly', 
    'Test Rate - Cumulative':'TestRateCumulative',
    'Percent Tested Positive - Weekly':'PercentTestedPositiveWeekly',
    'Percent Tested Positive - Cumulative':'PercentTestedPositiveCumulative', 
    'Deaths - Weekly':'DeathsWeekly',
    'Deaths - Cumulative':'DeathsCumulative',
    'Death Rate - Weekly':'DeathRateWeekly', 
    'Death Rate - Cumulative':'DeathRateCumulative',
    'Population':'Population', 
    'Row ID':'RowID',
    'x':'x', 
    'y':'y'      
    })  
    covid19_df=Covid_19_MultiIdx_weekZip.drop(['ZIP Code Location'], axis=1).copy()    
    Covid_19_MultiIdx_weekZip["geometry"]=Covid_19_MultiIdx_weekZip.apply(lambda row,x:shapely.wkt.loads(row[x]) if isinstance(row[x],str) else float('nan'),x='ZIP Code Location',axis=1)
    # print(Covid_19_MultiIdx_weekZip.xs(10,level=0,drop_level=False))
    # print(Covid_19_MultiIdx_weekZip.xs('60602',level=1,drop_level=False))
    
    Covid_19_MultiIdx_weekZip["x"]=Covid_19_MultiIdx_weekZip.apply(getPointCoords, geom='geometry', coord_type='x', axis=1)
    Covid_19_MultiIdx_weekZip["y"]=Covid_19_MultiIdx_weekZip.apply(getPointCoords, geom='geometry', coord_type='y', axis=1)
    # print(Covid_19_MultiIdx_weekZip.head())
    # print(Covid_19_MultiIdx_weekZip.dtypes)
    # print(Covid_19_MultiIdx_weekZip.index.values)
    # print(Covid_19_MultiIdx_weekZip.index)
    zip_codes= gpd.read_file(dataFpDic["zip_codes"])
    # print("-"*50)
    # print(zip_codes.columns)
    # print(zip_codes.head())
    covid19_df_joinColumn=covid19_df.reset_index(level=['Week Number','ZIP Code']).rename(columns={'ZIP Code':'zip'})
    covid19_zip=zip_codes.merge(covid19_df_joinColumn,on='zip')
        
    return Covid_19_MultiIdx_weekZip,covid19_df,covid19_zip,covid19_df_byZip

#显示vector数据