The issue you are having is due to axis order changes in PROJ 6+: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6
One solution is to use always_xy:
from shapely.geometry import LineString
from shapely.ops import transform
from functools import partial
roadsegment = LineString([(5.318945751388698, 50....
A simple approach is to use shapely.ops.transform, because that function is smart about using the return values to determine whether to create 2D or 3D shapes. This approach will however allocate new memory, so while efficient in code, it is not optimal for memory usage.
p = Point(2, 3) # no z
pz = shapely.ops.transform(lambda x, y: (x,...
With the help of @Taras I could solve my efficiency problem.
As the first try, I thought maybe the tqdm library is causing this delay. So I modified my code to
within_points = [points.geometry[i].within(map) for i in range(points.geometry.count())]
As I wasn't using tqdm to know the remaining time, I waited about one day to see if the process terminates ...
I got back a very helpful response from my County (Wake, NC) GIS help desk. Brandon guided me to the solution: "As to your other point, the multi-part polygon issue could be solved by dissolving the shapes (Dissolve in ArcGIS or Dissolve in GeoPandas). Dissolve will combine disparate polygons together based on an attribute field – in this case, the PIN_NUM ...
You could try rioxarray. An example of clipping/masking raster is here.
rds = rioxarray.open_rasterio('raster.asc')
geometries = [
geometry_crs = "???" # not sure what projection the geometries are in
rds_clipped = rds.rio.clip(...
@jelleve, your answer is correct and in line with https://pyproj4.github.io/pyproj/stable/advanced_examples.html#optimize-transformations.
Additionally, you may be interested in geopandas. It will enable you to read/write these geometries and reproject them using the same method you showed above using one line of code:
gdf = geopandas.read_file("input....
After some experimenting, it seems that I was just plain wrong. My initial code looked as follows:
transformer = pyproj.Transformer.from_proj(pyproj.Proj(init='epsg:27563'), pyproj.Proj(init='epsg:4326'))
projected = shapely.ops.transform(transformer.transform, polygon)
for polygon in polygons:
This is actually pretty easy.
Have a look at https://shapely.readthedocs.io/en/stable/manual.html#object.interpolate
distance = 0 ## starting at length 0
add_distance = 10 ## eg interpoalte every 10 meter
while distance < line.length:
new_point = line.interpolate(distance)
distance += add_distance ## add more
You will have to ...