Python shapely.geometry.mapping() 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 test_make_geocube__group_by_no_measurements(input_geodata, tmpdir):
    out_grid = make_geocube(
        vector_data=input_geodata,
        output_crs=TEST_GARS_PROJ,
        geom=json.dumps(mapping(TEST_GARS_POLY)),
        group_by="hzdept_r",
        resolution=(-10, 10),
        fill=-9999.0,
    )

    # test writing to netCDF
    out_grid.to_netcdf(
        str(tmpdir.mkdir("make_geocube_soil").join("soil_grid_group.nc"))
    )

    # test output data
    with xarray.open_dataset(
        os.path.join(TEST_COMPARE_DATA_DIR, "soil_grid_group.nc"), mask_and_scale=False
    ) as xdc:
        xarray.testing.assert_allclose(out_grid, xdc)

    tmpdir.remove() 

def __getitem__(self, geometry):
        if isinstance(geometry, BaseGeometry) or getattr(geometry, "__geo_interface__", None) is not None:
            if self._tms_meta._bounds is None:
                return self.aoi(geojson=mapping(geometry), from_proj=self.proj)
            image = GeoDaskImage.__getitem__(self, geometry)
            image._tms_meta = self._tms_meta
            return image
        else:
            result = super(TmsImage, self).__getitem__(geometry)
            image = super(TmsImage, self.__class__).__new__(self.__class__, result)
            if all([isinstance(e, slice) for e in geometry]) and len(geometry) == len(self.shape):
                xmin, ymin, xmax, ymax = geometry[2].start, geometry[1].start, geometry[2].stop, geometry[1].stop
                xmin = 0 if xmin is None else xmin
                ymin = 0 if ymin is None else ymin
                xmax = self.shape[2] if xmax is None else xmax
                ymax = self.shape[1] if ymax is None else ymax

                g = ops.transform(self.__geo_transform__.fwd, box(xmin, ymin, xmax, ymax))
                image.__geo_interface__ = mapping(g)
                image.__geo_transform__ = self.__geo_transform__ + (xmin, ymin)
            else:
                image.__geo_interface__ = self.__geo_interface__
                image.__geo_transform__ = self.__geo_transform__
            image._tms_meta = self._tms_meta
            return image 

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'}}) 

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 test_make_geocube__no_measurements(input_geodata, tmpdir):
    out_grid = make_geocube(
        vector_data=input_geodata,
        output_crs=TEST_GARS_PROJ,
        geom=json.dumps(mapping(TEST_GARS_POLY)),
        resolution=(-10, 10),
        fill=-9999.0,
    )

    # test writing to netCDF
    out_grid.to_netcdf(str(tmpdir.mkdir("make_geocube_soil").join("soil_grid_flat.nc")))

    # test output data
    with xarray.open_dataset(
        os.path.join(TEST_COMPARE_DATA_DIR, "soil_grid_flat.nc"), mask_and_scale=False
    ) as xdc:
        xarray.testing.assert_allclose(out_grid, xdc)

    tmpdir.remove() 

def __getattr__(self, name) -> Dict[str, Any]:
        if not name.startswith("osm_"):
            raise AttributeError(
                f"'{self.__class__.__name__}' has no attribute '{name}'"
            )
        else:
            values = self.osm_request().ways.values()
            return {
                "type": "FeatureCollection",
                "features": [
                    {
                        "type": "Feature",
                        "geometry": mapping(elt[1]),
                        "properties": elt[0]["tags"],
                    }
                    for elt in values
                    if elt[0]["tags"]["aeroway"] == name[4:]
                ],
            } 

def buffer_shape(infile, outfile, buffer=None):

    with collection(infile, "r") as in_shape:
        # schema = in_shape.schema.copy()
        schema = {'geometry': 'Polygon', 'properties': {'id': 'int'}}
        crs = in_shape.crs
        with collection(
                outfile, "w", "ESRI Shapefile", schema, crs=crs) as output:

            for i, point in enumerate(in_shape):
                output.write({
                    'properties': {
                        'id': i
                    },
                    'geometry': mapping(
                        shape(point['geometry']).buffer(buffer))
                }) 

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 

def make_record(properties: Dict[str, type], load_num: str, tree_row: Dict[str, str], point: Dict[str, str]):
    """
    shp 파일에 입력할 수목 정보 생성
    :return:
    """
    # 레코드 기본 구조
    record = {'properties': {"탐방로": load_num}, 'geometry': None}
    # 수목 정보 입력
    for key, value in zip(tree_row.keys(), tree_row.values()):
        atr_type = properties[key]
        # 속성 타입이 int인데 속성값이 ''일 경우 0으로 입력
        record['properties'][key] = atr_type(value) if value or atr_type is str else 0

    # 위치정보 입력
    record['properties']['경도'] = point['경도']
    record['properties']['위도'] = point['위도']
    record['properties']['고도'] = point['고도']
    record['geometry'] = mapping(Point(point['경도'], point['위도']))
    return record 

def get_final_value(self, value, as_json=False):
        json_value = format_geojson(mapping(value))
        if value.is_empty:
            return json_value
        rounded_value = shape(json_value)

        shapely_logger.setLevel('ERROR')
        if rounded_value.is_valid:
            return json_value if as_json else rounded_value
        shapely_logger.setLevel('INFO')

        rounded_value = rounded_value.buffer(0)
        if not rounded_value.is_empty:
            value = rounded_value
        else:
            logging.debug('Fixing rounded geometry failed, saving it to the database without rounding.')

        return format_geojson(mapping(value), rounded=False) if as_json else value 

def projectShapes(features, toCRS):
    import pyproj
    from functools import partial
    import fiona.crs as fcrs
    from shapely.geometry import shape, mapping
    from shapely.ops import transform as shpTrans

    project = partial(
        pyproj.transform,
        pyproj.Proj(fcrs.from_epsg(4326)),
        pyproj.Proj(toCRS))

    return list(
        {'geometry': mapping(
            shpTrans(
                project,
                shape(feat['geometry']))
        )} for feat in features
        ) 

def _load_CA_climate_zone_metadata(download_path):
    from shapely.geometry import shape, mapping

    ca_climate_zone_names = {
        "01": "Arcata",
        "02": "Santa Rosa",
        "03": "Oakland",
        "04": "San Jose-Reid",
        "05": "Santa Maria",
        "06": "Torrance",
        "07": "San Diego-Lindbergh",
        "08": "Fullerton",
        "09": "Burbank-Glendale",
        "10": "Riverside",
        "11": "Red Bluff",
        "12": "Sacramento",
        "13": "Fresno",
        "14": "Palmdale",
        "15": "Palm Spring-Intl",
        "16": "Blue Canyon",
    }
    geojson_path = os.path.join(
        download_path, "CA_Building_Standards_Climate_Zones.json"
    )
    with open(geojson_path, "r") as f:
        geojson = json.load(f)

    metadata = {}
    for feature in geojson["features"]:
        zone = "{:02d}".format(int(feature["properties"]["Zone"]))
        geometry = feature["geometry"]
        polygon = shape(geometry)
        metadata[zone] = {
            "ca_climate_zone": "CA_{}".format(zone),
            "name": ca_climate_zone_names[zone],
            "polygon": polygon,
            "geometry": to_geojson(polygon),
        }
    return metadata 

def export_tree_crowns(self, crowntype='crown_poly_smooth'):
        """ Convert tree crown raster to georeferenced polygon shapefile

        Parameters
        ----------
        crowntype :   str, optional
                      choose whether the raster of smoothed version should be
                      exported: `crown_poly_smooth` or `crown_poly_raster`
        """
        outfile = self.outpath / f'tree_{crowntype}.shp'
        outfile.parent.mkdir(parents=True, exist_ok=True)

        if outfile.exists():
            outfile.unlink()

        schema = {
            'geometry': 'Polygon',
            'properties': {'DN': 'int', 'TTH': 'float', 'TCH': 'float'}
        }
        with fiona.collection(
            str(outfile), 'w', 'ESRI Shapefile',
            schema, crs=self.srs
        ) as output:
            for tidx in range(len(self.trees)):
                feat = {}
                tree = self.trees.iloc[tidx]
                feat['geometry'] = mapping(tree[crowntype])
                feat['properties'] = {
                    'DN': tidx,
                    'TTH': float(tree.top_height),
                    'TCH': float(tree.top_cor_height)
                }
                output.write(feat) 

def export_tree_locations(self, loc='top'):
        """ Convert tree top raster indices to georeferenced 3D point shapefile

        Parameters
        ----------
        loc :     str, optional
                  tree seed position: `top` or `top_cor`
        """
        outfile = self.outpath / f'tree_location_{loc}.shp'
        outfile.parent.mkdir(parents=True, exist_ok=True)

        if outfile.exists():
            outfile.unlink()

        schema = {
            'geometry': '3D Point',
            'properties': {'DN': 'int', 'TH': 'float', 'COR': 'int'}
        }
        with fiona.collection(
            str(outfile), 'w', 'ESRI Shapefile', schema, crs=self.srs
        ) as output:
            for tidx in range(len(self.trees)):
                feat = {}
                tree = self.trees.iloc[tidx]
                feat['geometry'] = mapping(
                    Point(tree[loc].x, tree[loc].y, tree[f'{loc}_elevation'])
                )
                feat['properties'] = {'DN': tidx,
                                      'TH': float(tree[f'{loc}_height']),
                                      'COR': int(tree.tt_corrected)}
                output.write(feat) 

def get_feature_centroid(self, feature_geometry):
        """Gets the centroid from a feature's geometry"""
        if feature_geometry['type'].lower() == 'point':
            return feature_geometry['coordinates']
        else:
            geom = shape(feature_geometry)
            g_centroid = geom.centroid
            centroid = mapping(g_centroid)
            return centroid['coordinates'] 

def to_geojson(polygon):
    import simplejson
    from shapely.geometry import mapping

    return simplejson.dumps(pretty_floats(mapping(polygon)), separators=(",", ":")) 

def build_feature(tile_id, be_precise=None):
    """
    Given an SRTM tile ID, a list of (decoded) GeoJSON Feature object corresponding to the WGS84 longitude-latitude boundary of the tile.
    Use the same ``be_precise`` keyword as in :func:`get_bounds`.
    """
    return {
        'type': 'Feature',
        'properties': {'tile_id': tile_id},
        'geometry': mapping(build_polygon(tile_id, be_precise))
        } 

def _serialize(self, geometry=True, **kwargs):
        result = super()._serialize(geometry=geometry, **kwargs)
        result['width'] = float(str(self.width))
        result['height'] = float(str(self.height))
        result['altitude'] = float(str(self.altitude))
        result['color'] = self.color
        if geometry:
            result['buffered_geometry'] = format_geojson(mapping(self.buffered_geometry))
        return result 

def getGJSONinfo(geoJSONinfo):
    """
    Loads a lattice of GeoJSON, bounds, and creates a list mapping an on-the-fly UID w/ the actual index value.
    """
    features = list(i for i in filterBadJSON(geoJSONinfo))
    UIDs = list(feat['properties']['qt'] for feat in features)

    featDimensions = int(np.sqrt(len(features)/2.0))
    
    return features, UIDs, featDimensions 

def smart_mapping(geometry):
    if hasattr(geometry, 'wrapped_geojson'):
        return geometry.wrapped_geojson
    return shapely_mapping(geometry) 

def as_geojson(self):
        """ Returns a GeoJSON Feature object for this Destination. """
        return {
            "type": "Feature",
            "properties": {
                "name": self.name,
                "address": self.address,
            },
            "geometry": mapping(to_shape(self.point))
        } 

def make_centroid(self, geojson_geometry_str, geom_type='Polygon'):
        """ converts a geojson geojson_geometry_string
            OR a coordinate string into
            a centroid
        """
        centroid = False
        geojson_geom = False
        try:
            json_obj = json.loads(geojson_geometry_str)
        except:
            json_obj = False
        if isinstance(json_obj, list):
            geojson_geom = {'type': geom_type}
            geojson_geom['coordinates'] = json_obj
        elif isinstance(json_obj, dict):
            if 'type' in json_obj \
               and 'coordinates' in json_obj:
                geojson_geom = json_obj
            elif 'coordinates' in json_obj:
                geojson_geom = json_obj
                geojson_geom['type'] = geom_type
        else:
            geojson_geom = False
        if isinstance(geojson_geom, dict):
            geom = shape(geojson_geom)
            g_centroid = geom.centroid
            centroid = mapping(g_centroid)
        return centroid 

def insert(self,  # type: ignore
               keys: Union[Sequence[str], Mapping[str, str]],
               filepath: str, *,
               metadata: Mapping[str, Any] = None,
               skip_metadata: bool = False,
               override_path: str = None) -> None:
        """Insert a raster file into the database.

        Arguments:

            keys: Keys identifying the new dataset. Can either be given as a sequence of key
                values, or as a mapping ``{key_name: key_value}``.
            filepath: Path to the GDAL-readable raster file.
            metadata: If not given (default), call :meth:`compute_metadata` with default arguments
                to compute raster metadata. Otherwise, use the given values. This can be used to
                decouple metadata computation from insertion, or to use the optional arguments
                of :meth:`compute_metadata`.
            skip_metadata: Do not compute any raster metadata (will be computed during the first
                request instead). Use sparingly; this option has a detrimental result on the end
                user experience and might lead to surprising results. Has no effect if ``metadata``
                is given.
            override_path: Override the path to the raster file in the database. Use this option if
                you intend to copy the data somewhere else after insertion (e.g. when moving files
                to a cloud storage later on).

        """
        pass

    # specify signature and docstring for get_datasets 

def flight_kepler(flight: "Flight") -> Dict[str, Any]:
    return dict(
        geometry=mapping(flight.shape),
        properties={
            "icao24": flight.icao24,
            "callsign": flight.callsign,
            "registration": flight.registration,
            "start": f"{flight.start:%Y-%m-%d %H:%M:%S}",
            "stop": f"{flight.stop:%Y-%m-%d %H:%M:%S}",
        },
        type="Feature",
    ) 

def geojson(self):
        return {
            "type": "FeatureCollection",
            "features": [
                {
                    "geometry": mapping(shape),
                    "properties": info["tags"],
                    "type": "Feature",
                }
                for info, shape in self.osm_request().ways.values()
                if info["tags"]["aeroway"] != "aerodrome"
            ],
        } 

def geojson(self):
        """Returns the GeoJSON representation of the shape as a Dict.
        The transformation is delegated to shapely ``mapping`` method.
        """
        return mapping(self.shape) 

def export_json(self) -> Dict[str, Any]:
        export: Dict[str, Any] = {"name": self.name, "type": self.type}
        shapes = []
        for p in self:
            shapes.append(
                {
                    "upper": p.upper,
                    "lower": p.lower,
                    "polygon": mapping(p.polygon),
                }
            )
        export["shapes"] = shapes
        return export 

def query(self, searchAreaWkt, query, count=100, ttl='10s', index=default_index):
        '''
        Perform a vector services query using the QUERY API
        (https://gbdxdocs.digitalglobe.com/docs/vs-query-list-vector-items-returns-default-fields)

        ElasticSearch spatial indexing has some slop in it and can return some features that are 
        near to but not overlapping the search geometry. If you need precise overlapping of the
        search API you will need to run a geometric check on each result.
        
        If the caller requests more than 1000 records and it's possible that it will take longer than
        the default TTL value to pull a single page of 1000 records into memory, it's possible to raise
        the TTL duration by setting the 'ttl' parameter to something higher than the default of 10 seconds.
        For example, to set the TTL to 30 seconds, use '30s'.  For one minute, use '1m'.

        Args:
            searchAreaWkt: WKT Polygon of area to search
            query: Elastic Search query
            count: Maximum number of results to return, default is 100
            ttl: Amount of time for each temporary vector page to exist

        Returns:
            List of vector results
    
        '''
        if count < 1000:
            # issue a single page query
            search_area_polygon = load_wkt(searchAreaWkt)
            geojson = json.dumps(mapping(search_area_polygon))

            params = {
                "q": query,
                "count": min(count,1000),
            }

            url = self.query_index_url % index if index else self.query_url
            r = self.gbdx_connection.post(url, data=geojson, params=params)
            r.raise_for_status()
            return r.json()

        else:
            return list(self.query_iteratively(searchAreaWkt, query, count, ttl, index)) 

def _key_dict_to_sequence(self, keys: Union[Mapping[str, Any], Sequence[Any]]) -> List[Any]:
        """Convert {key_name: key_value} to [key_value] with the correct key order."""
        try:
            keys_as_mapping = cast(Mapping[str, Any], keys)
            return [keys_as_mapping[key] for key in self.key_names]
        except TypeError:  # not a mapping
            return list(keys)
        except KeyError as exc:
            raise exceptions.InvalidKeyError('Encountered unknown key') from exc