1

I am trying to creating a point shapefile that lies within polygon(a feature layer). By creating point geometry objects and filtering the ones that fall within the polygon(feature layer), I created a list (list of point geometry objects as "pointList"). Using CopyFeature_management tool, I created a shapefile out of "pointList". But this doesnot have projection and the attributes has FID, Shape and ID(all 0).

How could I project this point shapefile, the same projection as that of polygon?

Here is what I have so far, (where the extent is calculated for polygon feature layer):

pointList =[]
m=0
     x = random.uniform(extent.XMin, extent.XMax)
     y = random.uniform(extent.YMin, extent.YMax)
     points = arcpy.Point(x, y)
     pointGeo = arcpy.PointGeometry(points)  #creates point geometry object
     polygon = arcpy.Polygon(polygon_array)  #polygon object
     pointWithin = pointGeo.within(polygon)  #gives boolean output. 
     if pointWithin == True:
         m = m + 1
         pointList.append(pointGeo)
arcpy.CopyFeatures_management(pointList, outpath)
  • You will have to abandon copyfeatures. Get the CRS from your polygon feature class with arcpy.Describe(polygon).spatialReference, create a new feature class with CreateFeatureclass_management (insert CRS here), add your intended fields (X_coord, Y_coord) and use an arcpy.da.InsertCursor on your new feature class with ['shape@XY','X_coord','Y_coord'] as fields then insert new rows with ((x,y),x,y) if pointWithin - the first element is a tuple of X and Y then values for the fields X_Coord and Y_Coord. Or you can use fields ['Shape@','X_coord','Y_coord'] and insert (points,x,y) on your cursor.. – Michael Stimson Dec 10 '18 at 4:50
  • Do I start arcpy.da.InsertCursor (for new point featureclass )inplace of line where I am appending the list after the True condition? Doing so is giving me syntax error for that particular line. Thanks much for the feedback @MichaelStimson. – nima Dec 10 '18 at 16:35
1

I would do it like this:

import os, arcpy
SR = arcpy.SpatialReference(EPSG_Code)     # set this to match your data
outpath = r'c:\your\path\to\shapefile.shp' # set this to match your data

# create an empty point feature class
arcpy.CreateFeatureclass_management(os.path.dirname(outpath),os.path.basename(outpath),"POINT",spatial_reference=SR)
# add the desired fields to store your coordinate info
arcpy.AddField_management(outpath,'X_Coord','DOUBLE')
arcpy.AddField_management(outpath,'Y_Coord','DOUBLE')

# if your polygon array doesn't change then set it here
# rather than rebuilding on each iteration
polygon = arcpy.Polygon(polygon_array)  #polygon object
polygon.spatial_reference = SR          # assuming the same spatial reference

with arcpy.da.InsertCursor(outpath,['SHAPE@XY','X_Coord','Y_Coord']) as ICur:
    x = random.uniform(extent.XMin, extent.XMax)
    y = random.uniform(extent.YMin, extent.YMax)
    point       = arcpy.Point(x, y)
    pointGeo    = arcpy.PointGeometry(points)  #creates point geometry object
    pointGeo.spatial_reference = SR            # set the spatial reference of this point
    if pointGeo.within(polygon):
        # we want to keep this one because it's within
        # the focus polygon, now insert this row
        ICur.insertRow( ((x,y),x,y) )
        # note the brackets are important, the first is for insertRow
        # the second is for the row tuple and the third is for the X,Y
        # tuple required by SHAPE@XY.

The projection is defined when you create the shapefile (or GDB feature class) using Create Feature Class. Then establish an insert cursor (arcpy.da) in a with block so that the objects and locks will be cleaned up automatically on exit; in the help it states

When using InsertCursor on a point feature class, creating a PointGeometry and setting it to the SHAPE@ token is a comparatively expensive operation. Instead, define the point feature using tokens such as SHAPE@XY, SHAPE@Z, and SHAPE@M for faster, more efficient access.

So although we're creating a point geometry and setting its spatial reference it will only be used to test for within, if speed isn't of interest then by all means declare your insert cursor as with arcpy.da.InsertCursor(outpath,['SHAPE@','X_Coord','Y_Coord']) as ICur: then your insert row would be ICur.insertRow( (pointGeo,x,y) ).

If your points are in a different spatial reference to your polygons you can use the PointGeometry.projectAs(different_spatial_reference_object) to project to match your polygons; this method is valid for PointGeometry objects and not Point objects.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.