New answers tagged

1

Solution is fairly simple, maybe it'll be useful for somebody. Thanks to @martinfleis and this answer. First - we have point: point = Point(-1.7063,55.4127) Next create a circle basing on point and with buffer - in here 1000 (meters!): local_azimuthal_projection = f"+proj=aeqd +R=6371000 +units=m +lat_0={point.y} +lon_0={point.x}" wgs84_to_aeqd ...


2

All geometric operations that GeoPandas performs depends on Shapely package. Shapely performs all operations in the x-y plane. That means Shapely doesn't care if a shape is located at equator or at a pole, or it has geographic or projected coordinate system. 0.25 means 0.25 squared degree (0.25 deg2) which is meaningless in terms of area calculation. If ...


1

You have two alternatives. On the one hand, you can create/read a table from a PostGIS query using ST_Makeline using something like: SELECT ST_MakeLine(the_geom_webmercator ORDER BY track_seg_id ASC) AS the_geom_webmercator, track_fid, row_number() OVER() AS cartodb_id FROM table_name GROUP BY track_fid On the other,...


3

I reported this as a bug here: https://github.com/pyproj4/pyproj/issues/549 In the meantime, you can keep using the deprecated "{'init':'epsg:4326'}" syntax a bit longer (and ignore the warning for now). Alternatively, you can use the following workaround to read the CRS directly from the .prj file (avoiding going through GDAL): gdf = geopandas.read_file('...


2

You can use this script: import utm import geopandas as gpd df = gpd.read_file("PATH/TO/Pipes.shp") def utm_to_latlon(coords, zone_number, zone_letter): easting = coords[0] northing = coords[1] return utm.to_latlon(easting, northing, zone_number, zone_letter) # Using nested list comprehension df ["lat_lon_tuple"] = [[utm_to_latlon(xy, 44, "N")...


2

Try using a Spatial Join with GeoPandas. If you wish you can subset the results to show only the points contained within your polygons using notnull(). This example performs a spatial join with 1,000,000 points and 10 polygons (i.e. A-J) in 42 seconds. import geopandas as gpd points_shp = '/path/to/points.shp' polys_shp = '/path/to/polygons.shp' points = ...


1

You just need to create the bounding box of each feature you are creating. Take a look at this thread there is a sample on how to achieve that. Find bounding box for multiple features using PyQGIS? You could achieve that by creating top right point (Xmax,Ymax) and (Xmin,Ymin) and return the bounding box of the collection of the two points. This is the ...


1

In case of GeoDataFrame, CRS in GeoPandas is stored on the level of GeoDataFrame, not individual GeoSeries (as of version 0.7.0, there is a discussion to change it). At this moment, I think that your solution of reprojecting GeoSeries and then assigning then to GeoDataFrame is the best solution, although admittedly not very elegant. Feel free to express your ...


1

You could save the classification array as a raster in memory (using gdal's MEM driver to create it) and then use gdal.Polygonize() function. Here is an example of how to use it. If you definitely do not want to create a raster file (not even in memory) another approach would be using the rasterio.features.shapes() function. What this function does is: ...


0

I have a solution which is requires only psycopg2 and shapely (in addition geopandas of course). It is generally bad practice to iterate through (Geo)DataFrame objects because it is slow, but for small ones, or for one-off tasks, it will still get the job done. Basically it works by dumping the geometry to WKB format in another column and then re-casts it ...


2

Yes both are correct. The difference is because one it taking into account the curvature of the earth while the other is not. See following images: and


3

In geopandas <= v0.6.3., gdf.crs returns a dictionary like {'init': 'epsg:EPSG_CODE'}. So, more appropriate way is to use tools of geopandas defined in geopandas.tools module. geom_srid_num = gpd.tools.crs.epsg_from_crs(gdf.crs) print(geom_srid_num) # OUT: 32616 -> int EDIT: As @snowman2 states in comment, epsg_from_crs is deprecated in geopandas v0....


0

One possible way is as follows. I am sure though, there are nicer ways! def seedtoleech(seed_layer,leech_layer,seed_column): # First, the index is reset. # Both indices will be transfered to the gpd.overlay() layer # to allow identifying overlaps. seed_layer = seed_layer.reset_index() leech_layer = leech_layer.reset_index() #...


0

Perhaps shapely's unary_union will help, where df is your GeoDataFrame: import pandas as pd from shapely.ops import unary_union hole_indx = pd.isna(df.raster_val) holes = unary_union(df.loc[hole_indx, 'geometry']) polys = unary_union(df.loc[~hole_indx, 'geometry']) polys_and_holes = polys.difference(holes)


3

You should try to preserve the WKT form if possible. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems It all depends on the version of geopandas. When version 0.7.0 comes out, you can do: boxes = geopandas.GeoDataFrame() with rasterio.open(raster_path) as dataset: boxes.crs = dataset.crs boxes.to_file(...


0

This solution worked for me. Just adjust the buffer amount to work with your units. Solution copied here: polygon.buffer(10, join_style=1).buffer(-10.0, join_style=1)


0

Following the answer given by @sbphd, this is what I coded. import geopandas as gpd # with columns "id", "latitude", "longitude" - 10k records df gdf = gpd.GeoDataFrame( df, geometry=gpd.points_from_xy( df["longitude"], df["latitude"], ), crs={"init":"EPSG:4326"}, ) # 10 records filtered_df filtered_gdf = gpd.GeoDataFrame(...


1

geocube is a new tool specifically designed for rasterizing geopandas data that wraps rasterio. It simplifies the process and eliminates the need for a template raster. https://github.com/corteva/geocube In the context of the example above: from geocube.api.core import make_geocube import geopandas counties = geopandas.read_file("zip://...


2

Not really GIS related, but you can use the astype method: geo_df['denominator'] = df[["basalareap","basalareas","basalaread"]].sum(axis=1).astype(str)


1

Create GeoPandas geodataframes: import geopandas as gpd import shapely df # your pandas dataframe with 10k records df_filt # your filtered dataframe # Create geometries from your lat-lons geom_list = [shapely.geometry.Point(lon,lat) for lon,lat in zip(df["longitude" ,df["latitude"])] # check the ordering of lon/lat # create geopandas geodataframe gdf = ...


3

You may be interested in geocube (https://github.com/corteva/geocube): Examples: https://corteva.github.io/geocube/stable/examples/examples.html from geocube.api.core import make_geocube gdf = geopandas.read_file(...) cube = make_geocube(vector_data=gdf, measurements=["column_name"], resolution=(16, -16)) out_grid["column_name"].rio.to_raster("...


4

You can first union all polygons with unary_union: single_multi_polygon = all_Rapa_Nui.unary_union This should now be a single MultiPolygon consisting of two polygons for the two islands. And then you can get the polygon parts of this MultiPolygon: polygons = single_multi_polygon.geoms


1

You should specify a common key using on parameter. I removed merged_ prefix for legibility. df = spatial_df.merge(tab_df, on='mukey', how='left') # df = tab_df.merge(spatial_df, on='mukey', how='right') gdf = gpd.GeoDataFrame(df) Sample spatial_df: col1 col2 mukey geometry 0 A 1.76 1 ... 1 B 0.40 2 ... 2 C 0.97 3 ... ...


6

I tried several coordinate systems for Florida and the various State Plane Florida East zone but the coordinates just didn't seem to fit. Google kept finding the address in Hempstead, NY, so I tried the Long Island zone, unit of US survey feet, and the results look good. Try EPSG::6539 for the NAD83 (2011) version or EPSG::2908 for NAD83 (HARN). For the ...


Top 50 recent answers are included