I'm trying to make a spatial verification. I have a geodataframe with 245 points, and I'm having a bad time trying to check if they're inside a polygon or not. I've been attempting this code, but I just won't work. As far as I understood, geodataframes are slightly different from dataframes when in a for loop. for x in event_log_c.geometry:     i = event_log_c.iloc[x].within(hexa)     print(i)  After running this, it brings me this error: TypeError: Cannot index by location index with a non-integer key  hexa is the polygon geodataframe, and event_log_c is the 245 points geodataframe. How can I make a spatial verification iterating a geodataframe? Are there any way of doing this?
122454
I'm attempting to do a spatial join much like the example here: . However, that approach seems really inefficient / slow. Even running this with a modest 250 points takes almost 2 minutes and it fails entirely on shapefiles with > 1,000 points. Is there a better approach? I'd like to do this entirely in Python without using ArcGIS, QGIS, etc.  I'd also be interested to know if it's possible to SUM attributes (i.e. population) of all the points that fall within a polygon and join that quantity to the polygon shapefile.       Here is the code I'm trying to convert. I get an error on line 9:  poly['properties']['score'] += point['properties']['score']   which says:     TypeError: unsupported operand type(s) for +=: 'NoneType' and 'float'.    If I replace the "+=" with "=" it runs fine but that doesn't sum the fields. I've also tried making these as integers but that fails as well.  with fiona.open(poly_shp, 'r') as n:    with fiona.open(point_shp,'r') as s:     outSchema = {'geometry': 'Polygon','properties':{'region':'str','score':'float'}}     with fiona.open (out_shp, 'w', 'ESRI Shapefile', outSchema, crs) as output:         for point in s:             for poly in n:                 if shape(point['geometry']).within(shape(poly['geometry'])):                       poly['properties']['score']) += point['properties']['score'])                     output.write({                         'properties':{                             'region':poly['properties']['NAME'],                             'score':poly['properties']['score']},                         'geometry':poly['geometry']})
112601
I've just posted this at StackOverflow and someone suggested I try here. Original:   I'm trying to represent a rectangular area which crosses 180 degrees longitude. For more background see   Here's my test case:  from django.contrib.gis.geos import Polygon, MultiPolygon from my_project.my_app.models import Photo  a = Polygon.from_bbox((30, -80, 180, 80))  # the part to the east of 180 b = Polygon.from_bbox((-180, -80, -160, 80))  # a part to the west of 180 c = Polygon.from_bbox((-180, -80, -100, 80))  # a larger part to the west of 180  ok   = MultiPolygon(a,b) ok2  = MultiPolygon(c) boom = MultiPolygon(a,c)  # This works Photo.objects.filter(location__coveredby=ok)[:1] # This also works so c is ok Photo.objects.filter(location__coveredby=ok2)[:1] # This gives "BOOM! Could not generate outside point!" Photo.objects.filter(location__coveredby=boom)[:1]  # splitting c doesn't help c1 = Polygon.from_bbox((-180, -80, -140, 80)) c2 = Polygon.from_bbox((-140, -80, -100, 80)) test = MultiPolygon(a,c1,c2) Photo.objects.filter(location__coveredby=test)[:1] # BOOM! Could not generate outside point!   By changing the numbers I can make this error come and go. (-180, -80, x, 80) works where x &lt;= -140 for example. For every number there is a threshold like this but I can't find a pattern. For boxes with the same area, some will work and others not. For boxes with the same width some will work and others not.   I can look at the SQL being generated but the areas are represented in binary (EWKB) and I'm not sure how to read it.  Can anyone explain this?
115683
Let me describe the data setup. A line has been provided as an array of coordinate pairs, one for each vertex. This array is variable in size, anywhere from 2 to 200 coordinate pairs per line. What I am trying to wrap my mind around is the best way to create these lines. Arcpy and ArcObjects are both viable options, but ArcObjects is preferred due to its speed (generally).   I have brainstormed a few methods, but they fall short for some reason or another:  1) I could use the  with arcpy to create polylines and then Copy Features, but the resulting geometry would not have any of the attribution necessary to identify said points.  2) I could use , but this only works with one pair of point, and I require something much more robust.  The long and short is that I require a method that can create a polyline from arrays and maintain a schema for data that I have created in a table. One idea I have is to create said table and schema, and then somehow create the polyline and copy its shape property to the output table, porting over the ID and whatever other attribution is needed. How can I go about creating the polyline's shape property without resorting to an intermediate shapefile? Or, if there is a better method that comes to mind, please share it with me.  Update This code seems to work!  IFeatureClass fClass = gpUtil.OpenFeatureClassFromString(path); IFeatureCursor cursor = fClass.Insert(true); IFeatureBuffer buffer = fClass.CreateFeatureBuffer();  &lt;snip...&gt;  string[] nodearray = a.Split(';'); foreach (string i in nodearray) {     query.WhereClause = "\"PNT_ID\" = '" + i+"'";     Console.WriteLine(query.WhereClause.ToString());     ICursor search = table.Search(query, false);     IRow row = search.NextRow();     float x = float.Parse(row.get_Value(2).ToString());     float y = float.Parse(row.get_Value(3).ToString());     int id = int.Parse(row.get_Value(1).ToString());     point.X = x;     point.Y = y;     point.ID = id;     point.SpatialReference = sRef;     pCollection.AddPoint(point, ref Missing, ref Missing);  } IGeometry5 polygeo = (IGeometry5)pCollection; buffer.set_Value(1, polygeo); buffer.set_Value(2, wayid); cursor.InsertFeature(buffer);
52459
I've a Google map with bunch of polygons on it.  Here is a problem I'm interested in: Given a lat,lng point, what's the best way to determine all the polygons that this point lies in?  The obvious way is to run a "point in polygon" algorithm iteratively for each polygon, but I was wondering if there an efficient algorithm to answer such queries especially if you have thousands of polygons.
607
I was asked in a recent training course if QGIS could automatically calculate the next/previous and above/below page numbers for a map book created using the atlas generator. I managed to work out a fairly reasonable label expression for a regular grid if you know the width and height of the grid.  But we then started to think of realistic examples where we don't want to draw pages that don't contain our district of interest, such as this one of my home county:    So this afternoon I had a play at a python script to work out the 4 neighbours I was interested in for each grid cell and added those values to my grid (this is heavily based on Ujaval Gandhi's ):  for f in feature_dict.values():     print 'Working on %s' % f[_NAME_FIELD]     geom = f.geometry()     # Find all features that intersect the bounding box of the current feature.     # We use spatial index to find the features intersecting the bounding box     # of the current feature. This will narrow down the features that we need     # to check neighboring features.     intersecting_ids = index.intersects(geom.boundingBox())     # Initalize neighbors list and sum     neighbors = []     neighbors_sum = 0     for intersecting_id in intersecting_ids:         # Look up the feature from the dictionary         intersecting_f = feature_dict[intersecting_id]         int_geom = intersecting_f.geometry()         centroid = geom.centroid()         height = geom.boundingBox().height()         width = geom.boundingBox().width()         # For our purpose we consider a feature as 'neighbor' if it touches or         # intersects a feature. We use the 'disjoint' predicate to satisfy         # these conditions. So if a feature is not disjoint, it is a neighbor.         if (f != intersecting_f and             not int_geom.disjoint(geom)):             above_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),                centroid.asPoint().y()+height))             below_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),                centroid.asPoint().y()-height))             left_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()-width,                centroid.asPoint().y()))             right_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()+width,                centroid.asPoint().y()))             above = int_geom.contains(above_point)                below = int_geom.contains(below_point)                left = int_geom.contains(left_point)             right = int_geom.contains(right_point)             if above:                 print "setting %d as above %d"%(intersecting_f['id'],f['id'])                 f['above']=intersecting_f['id']              if below:                 print "setting %d as below %d"%(intersecting_f['id'],f['id'])                 f['below']=intersecting_f['id']              if left:                 print "setting %d as left of %d"%(intersecting_f['id'],f['id'])                 f['left']=intersecting_f['id']              if right:                 print "setting %d as right of %d"%(intersecting_f['id'],f['id'])                 f['right']=intersecting_f['id']      # Update the layer with new attribute values.     layer.updateFeature(f)  layer.commitChanges()   This works just fine.     But to be honest the whole creating a test point to the North and then testing all the possible neighbours seems wrong. However after an afternoon of wracking my brain I can't think of a better way to determine what is a particular grid cell's northern neighbour?  Ideally I'd like something simple enough to put in a print composer text box, but I suspect that's too much to ask for.
43815
So I found this python script on another post here and was trying to adjust it for my needs.  I'm a very novice python user so I'm struggling with how to modify the script.  I have a Feature Class stored in a feature Dataset that I want to search a field for duplicate values and populate a new field with Y for duplicate or N for none.  The below script is what I found that looks like it will work once I find a way to drill down into my file geodatabase.          from arcpy import *          inShapefile = pointsShapefile          checkField = "xyCombine"         updateField = "dplicate"         #List of values found once         occursOnce = []        #list of values found twice         occursTwice = []       cursor = da.SearchCursor (inShapefile, [checkField])           for row in cursor:           #Check value is not null        if row[0]:     #If not already found to occur twice, proceed     if not row[0] in occursTwice:         #If hasn't occured once yet         if not row[0] in occursOnce:             #Add to occurs once list             occursOnce.append (row[0])         #If value has already been found once         else:             #Add to occurs twice list (duplicates)             occursTwice.append (row[0])     del cursor       cursor = da.UpdateCursor (inShapefile, [checkField, updateField])      for row in cursor:       #Check value is not null         if row[0]:            #check if value in occursTwice list (i.e. is duplicate)     if row[0] in occursTwice:         row[1] = "Y"     else:         row[1] = "N"     cursor.updateRow(row)    del cursor
115738
I am trying to import an Excel table with around 10,000 XY points into ArcMap 10.2. I am using State Plane NAD 1983 Alabama West (US Feet) as the Coordinate System in my map data frame. I then add the excel table to my Table of Contents in the map then select Display XY Data. Next, I choose the same coordinate system as my data frame then select OK. Once it creates the Events layer for my points they are not in the correct location.   Should I switch everything in my map to WGS84 or is there a solution for my points to be added onto my current map with the state plane coordinate system?
76274
I'm exploring how we can use NetTopologySuite to address a geoprocessing task. based on some code sample I found, I was able to piece together a sample application that reads a shapefile and checks if a point falls within a polygon. The issue is, my point is in WGS84 but the shapefile is in EPSG:2810. I need to transform the input point to the CRS of the shapefile.  How do I convert a point from WGS to EPSG:2810 using NetTopologySuite 1.13.3?  Here's some of the code I have so far:  string shpFilename = @"C:\...\polygons.shp"; GeometryFactory factory = new GeometryFactory(); var coord = new Coordinate(-90.00, 45.00); var point = factory.CreatePoint(coord);   I should also note that when I manually convert from WGS84 to the projected system and use those converted values, my code works as designed. This confirms this is a matter of transforming the coordinates before attempting to determine if they fall within a polygon.
27958
Is there a simple way / ready-to-go function to calculate the azimuth of 2D-polygons in a GeoPandas DataFrame? My goal is to eventually compute azimuth angles of building footprints. Most of the functions I have found explicitly calculate the azimuth for a line, which is not what I am looking for. I have been thinking along the following lines to calculate the azimuth of a polygon using its normal vector (2 dimensions). Is my approach going into the right directions? I am a bit unsure how to adjust for the quadrants, also my results seem a bit odd. azimuth = 90 - math.degrees(math.atan2(normal[1], normal[0])) if azimuth &gt;= 360.0:     azimuth -= 360.0 elif azimuth &lt; 0.0:     azimuth += 360.0
21645
the point represents rates like 10%, etc. (this was also pretty clear from 2 years ago) -- been trying to do this for 2 yrs...  the points dont represent numbers of events since they are rates ()    but the points have geo ids so they are attached to location, not distance (whatever distance is to imply) -- update: it represents z-value but we're just dealing with rates
69554
I am trying to get coordinates from a mouse click in QGIS 3 and python 3 using the plugin builder.  I have tried using the response , and i can manage to get the coordinates printed to the python console in QGIS but at the same time i get the message that "'TypeError: 'QgsMapToolEmitPoint' object is not iterable" or "'QgsMapToolEmitPoint' object does not support indexing"   I can't remember i ever seen a python program causing an error but still manage to compute and output the wanted result. The "not iterable" error traces back to the line with the for loop "for point in ... " the "does not support indexing" error is caused if i do e.g. x = pointTool[0]  Here is a version of the code:  from qgis.gui import QgsMapToolEmitPoint   def display_point(pointTool):     #print(pointTool)     cood = []     for point in pointTool:         cood.append(point)     print(cood)  # a reference to our map canvas canvas = iface.mapCanvas()  # this QGIS tool emits as QgsPoint after each click on the map canvas pointTool = QgsMapToolEmitPoint(canvas)  pointTool.canvasClicked.connect(display_point)  canvas.setMapTool(pointTool)  display_point(pointTool)
11680
My question is about extracting values from rasters at the location of points. With the function extract this is very easy, and the function gives me a dataframe with the values of all the variables in the points. I want to have in that dataframe the coordinates of each point. How can I make that happen? Is it possible to say to R that when extracting the values from the raster also add the columns of the location point.  This is my procedure:  presencias=read.table("c:/SDM_R/presencias/P_lentiscus_pres.csv",header=TRUE,sep=";")  lista_variables &lt;-list.files(path="Variables_modelizacion/solo_ascii",pattern='*.asc',full.names=TRUE) variables &lt;- stack(lista_variables)  variables_presencia&lt;-extract(variables, presencias)   the results are something like these:  &gt; bio1  bio12  bio18  bio2  bio4  &gt; 90    875    165    95    4886 &gt; 115   1085   158    83    4075 &gt; 135   1153   153    67    3402 &gt; 85    1026   137    99    5203 &gt; 96    667    128    108   5823 &gt; 98    531    109    113   6305 &gt; 132   450    63     123   6598 &gt; 132   569    104    106   5963 &gt; 95    814    196    98    5571 &gt; 146   474    39     114   6603   But I want two more columns with the coordinate data (but could be others columns from my csv extracting location table).  Thank so much.
14395
I am attempting to iterate over all the feature classes in a folder and I'm getting the error 000368. The other solutions in on SE were not helpful to me unfortunately, so I'm posting a separate question. search_grids is a list of point feature classes and searches is a another point feature class.   The points never actually intersect with each other, but I'm using the Select by location tool to select points within a specified distance of each other. This isn't the reason for the error though as I was able to get this to work in the python window in ArcMap.   My code is the following  arcpy.env.workspace = r"C:\workspace" search_grids = arcpy.ListFeatureClasses() searches = r"C:\workspace\search" for grid in search_grids:     arcpy.management.SelectLayerByLocation(grid,"INTERSECT",searches,"97.5 meters","NEW_SELECTION","INVERT")     arcpy.management.SelectLayerByLocation(grid,"INTERSECT",searches,"2.5 meters","ADD_TO_SELECTION","NOT_INVERT")   Why am I getting this error?       arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.   ERROR 000368: Invalid input data.   Failed to execute (SelectLayerByLocation)
85729
Is it possible (without writing Python code - sorry I'm not a developer) to create a point layer from a polygon layer but ensure the point is within the polygon? The ftools centroid tool is cool but some centroids are created outside the polygon (i.e. irregular polygon). In ArcGIS there is a "Feature to Point" tool similar to the ftools 'polygon centroids' tool but has a checkbox to force the point to be within the polygon.
78046
My eventual goal is to automate kernel density calculations for a large corpus of point data; however, if all points are identical (for example, if there are 7 points all stacked on top of one another because they have the same long/lat values), the  in arcpy will not work. It will result in the following error:  "Error 010246: All input points are identical A valid window must be specified"   Since this will be an automated process, I need a way to check the points in the current input to be sure they aren't all the same using python/arcpy. Is there a way to do that?   EDIT: I specified a grid size and an extent, so a lot of the errors are fixed; however, identical points still give me the "valid window must be specified." I have tried to use the FindIdentical function to look for identical points, but this returns a DBF file, that I cannot read in the script. Is there some way to read a DBF easily in python, or can I find if all points are identical in another way?
107585
I have an Excel file with the coordinates (Latitude/Longitude) for an area. The file contains the coordinate of two corners in the latitude/longitude format.  How can we do to use Python to iterate the process for all coordinates and to export into XY coordinates?  Example of my file:  N°1:  Min_Long   Max_Long  Min_Lat   Max_Lat      1.5        2.833      6.25         7.33  N°2  -13        5          10          25   ... ...  Or else is there a tool in ArcGIS similar to "Add X,Y Data" in ArcMap, that can be used to automatically add this lat/long format?
54545
I'm using PostGIS 1.5.2, with geometries in SRID:900913. I need create a circles using a list of points as centre, with a radius of 600 kilometers. I'm using this query:  INSERT INTO circles (geom) (    SELECT ST_Buffer(point, 600000, 'quad_segs=8')     FROM points );   But circles created have not 600 kilometers of radius (radius is near this length, but not exactly).  Are there other methods to create circles in PostGIS?  NOTE: Information represented is from Spain. Correct projection is 4326, but client use Google rasters, so I am storing data in 900913 to avoid reprojections and increase performance.