1

I am writing a script which scrapes USGS earthquake data in GeoJSON format. I then try to create a shapefile using Arcpy copy features tool and I cannot seem to get it right

import urllib2
import json
import arcpy
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.geojson"
weburl = urllib2.urlopen(url)
if weburl.getcode() == 200:
    data = json.loads(weburl.read())
earthquakes = {}
for i in data["features"]:
    mag, place = i["properties"]["mag"],i["properties"]["place"]
    x,y = i["geometry"]["coordinates"][0],i["geometry"]["coordinates"][1]
    earthquakes[place] = [mag,x,y]
pointlist = []
point = arcpy.Point()
for k,v in earthquakes.items():
    point.X = v[1]
    point.Y = v[2]
    magn = v[0]
    place = k
    pointGeometry = arcpy.PointGeometry(point)
    pointlist.append(place)
    pointlist.append(magn)
    pointlist.append(pointGeometry)
arcpy.CopyFeatures_management(pointlist, "C:\Users\rzagha\Desktop\earthquake.shp")

I keep getting this error

Traceback (most recent call last):
  File "<string>", line 254, in run_nodebug
  File "<module1>", line 24, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 2429, in CopyFeatures
    raise e
RuntimeError: Object: Error in executing tool

when I print out the pointlist:

[u'29km SSW of Leilani Estates, Hawaii', 2.65, <PointGeometry object at 0x7dd47f0L[0x7d6d530L]>, u'67km NNW of Jayapura, Indonesia', 4.7, <PointGeometry object at 0x7dd4198L[0x128f3490L]>, u'61km NNE of Isabela, Puerto Rico', 2.9, <PointGeometry object at 0x7dd4668L[0x128f3c88L]>, u'61km N of Isabela, Puerto Rico', 2.9, <PointGeometry object at 0x7dd4278L[0x5917968L]>, u'150km SSW of Ndoi Island, Fiji', 5, <PointGeometry object at 0x7dd4240L[0x5917af8L]>, u'19km NE of La Paz Centro, Nicaragua', 5.5, <PointGeometry object at 0x7dd49e8L[0x5917ad0L]>, u'11km ESE of East Foothills, California', 3.25, <PointGeometry object at 0x7dd4208L[0x5917b20L]>, u'134km SE of Pangai, Tonga', 5, <PointGeometry object at 0x7dd4470L[0x5917b70L]>, u'23km ESE of Three Rivers, Oregon', 2.8, <PointGeometry object at 0x7dd4a90L[0x12b42f80L]>, u'27km ENE of Winnemucca, Nevada', 2.6, <PointGeometry object at 0x7dd4320L[0x12b42e68L]>, u'66km E of Namie, Japan', 4.6, <PointGeometry object at 0x7dd41d0L[0x12b42d28L]>, u'23km ENE of Mooreland, Oklahoma', 2.5, <PointGeometry object at 0x7dd4518L[0x12b42f58L]>, u'44km WNW of Piru, Indonesia', 4.6, <PointGeometry object at 0x7dcbdd8L[0x12b42a80L]>, u'13km WNW of Pawnee, Oklahoma', 2.8, <PointGeometry object at 0x7dcba90L[0x12b42d00L]>, u'21km SSE of Ridgemark, California', 2.84, <PointGeometry object at 0x7dcb668L[0x12b42760L]>, u'20km ESE of Little Lake, CA', 3.06, <PointGeometry object at 0x7dcbb70L[0x12b42dc8L]>, u'78km E of Old Iliamna, Alaska', 2.5, <PointGeometry object at 0x7dcb978L[0x12b42ee0L]>, u'Southeast of Easter Island', 5.2, <PointGeometry object at 0x7dcbd30L[0x12b42670L]>, u'47km S of Kangding, China', 4.4, <PointGeometry object at 0x7dcbcf8L[0x12b42d78L]>, u'114km S of Raoul Island, New Zealand', 5.3, <PointGeometry object at 0x7dcb8d0L[0x12b42c88L]>, u'126km SSE of Lambasa, Fiji', 4.5, <PointGeometry object at 0x7dcb6a0L[0x12b42738L]>, u'192km ENE of Saipan, Northern Mariana Islands', 5, <PointGeometry object at 0x7dcbc50L[0x12b42710L]>, u'99km S of San Pedro de Atacama, Chile', 4.2, <PointGeometry object at 0x7dcbc88L[0x12b428f0L]>, u'139km E of Shikotan, Russia', 4.8, <PointGeometry object at 0x5944c88L[0x12b42e90L]>, u'61km SSW of Santiago Pinotepa Nacional, Mexico', 5.1, <PointGeometry object at 0x7d96978L[0x12b425f8L]>, u'10km ESE of Cuajinicuilapa, Mexico', 5.2, <PointGeometry object at 0x59189b0L[0x12b42418L]>, u'33km N of Sutton-Alpine, Alaska', 2.9, <PointGeometry object at 0x59187b8L[0x12b42c10L]>, u'54km ENE of Fritz Creek, Alaska', 3, <PointGeometry object at 0x5918710L[0x12b42fd0L]>, u'North of Ascension Island', 4.9, <PointGeometry object at 0x5918908L[0x12b42da0L]>, u'78km SE of Lata, Solomon Islands', 4.8, <PointGeometry object at 0xacff978L[0x12b48468L]>, u'20km SSE of Guanica, Puerto Rico', 2.5, <PointGeometry object at 0x12b44048L[0x12b48e90L]>]

I am using this ESRI page to guide me http://pro.arcgis.com/en/pro-app/arcpy/classes/pointgeometry.htm

Should I try and create a feature class then use an insert cursor to update that shapefile? what is the reason the copy features tool is not taking in my pointgeometry?

I have also tried removing the place and magn variables and used just the x,y coordinates and I get the same error


from @artwork21 answer: i created a feature class in a GDB and used an insert cursor to insert the rows

def GeoJSONToFC(gdb,USGSurl,fcname):
    import urllib2
    import json
    import arcpy
    arcpy.env.overwriteOutput = True
    geodb = gdb
    url = USGSurl
    weburl = urllib2.urlopen(url)
    if weburl.getcode() == 200:
        data = json.loads(weburl.read())
    earthquakes = []
    for i in data["features"]:
        mag, place = i["properties"]["mag"],i["properties"]["place"]
        x,y = float(i["geometry"]["coordinates"][0]),float(i["geometry"]["coordinates"][1])
        earthquakes.append([place,mag,[x,y]])
    earthq = arcpy.CreateFeatureclass_management(geodb,fcname,"POINT",'','',"ENABLED",4326)
    arcpy.AddField_management(earthq,"Place","TEXT")
    arcpy.AddField_management (earthq, "Magnitude", "SHORT")
    with arcpy.da.InsertCursor(earthq,["Place","Magnitude","SHAPE@XY"]) as cur:
        for v in earthquakes:
            place = v[0]
            mag = v[1]
            pointGeometry = arcpy.PointGeometry(arcpy.Point(*v[2]))
            row = (place,mag,pointGeometry)
            cur.insertRow(row)
            print row

gdb = "C:\Users\Ralph\Desktop\quakes.gdb"
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_month.geojson"
name = "MajorQuakes"
GeoJSONToFC(gdb,url,name)
  • the copy tool requires a Feature Layer, shapefile or feature class input format, which is why it is failing. You may want to look at arcpy.JSONToFeatures_conversion tool which will output a feature class or shapefile. ANd yes, pre-creating the feature class and using a cursor is another valid approach. – Ben S Nadler Sep 29 '16 at 15:31
  • @BenSNadler I am following the ESRI documentation on creating a Point Shapefile, once I scrape the values off the web page and store them into variables it should be no different than extracting coordinates from a csv or txt file – ziggy Sep 29 '16 at 15:33
3

Copy Features method takes in a Feature Layer, you are passing in an array. As you commented, you should create a empty feature class and use and insert cursor to create the feature and include the related field attribute values.

Copies features from the input feature class or layer to a new feature class. If the input is a layer which has a selection, only the selected features will be copied. If the input is a geodatabase feature class or shapefile, all features will be copied.

  • I surprisingly was able use copy features to just work on the XY coordinates but i was not able to get the place and magnitude fields in there. I took your advice and created a feature class, and used an insert cursor to update the field and it worked. il will udpate my question with the answer – ziggy Sep 29 '16 at 16:30
5

Just wanted to add that what you were trying to do does work with some small tweaks:

import urllib2
import json
import arcpy
url = "http://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.geojson"
weburl = urllib2.urlopen(url)
if weburl.getcode() == 200:
    data = json.loads(weburl.read())
earthquakes = []
for i in data["features"]:
    mag, place = i["properties"]["mag"],i["properties"]["place"]
    x,y = i["geometry"]["coordinates"][0],i["geometry"]["coordinates"][1]
    earthquakes.append([place,mag,[x,y]])
m = []
sr = arcpy.SpatialReference(4326)
for v in earthquakes:
    place = v[0]
    magn = v[1]
    m.append(arcpy.PointGeometry(arcpy.Point(*v[2]), sr))
    print place,magn,type(m)
arcpy.CopyFeatures_management(m, "C:\\Users\\rzagha\\Desktop\\earthquake.shp")

For copy features to work this way, you can pass in a list of arcpy.PointGeometry() objects. With this, you do need to provide a spatial reference (wkid 4326 is WGS1984).

However, as @artwork21 suggested, Creating a new feature class with all the schema predefined with an insert cursor would be the way to go. Then you could just map all the geojson field info in there at the same time.

Also, it is worth noting you can get a PointGeometry object by using the arcpy.AsShape() function:

pt = arcpy.AsShape(i["geometry"], False)
pt.spatial_reference = arcpy.SpatialReference(4326) # set sr to WGS 1984

This is the same as:

pt = arcpy.AsShape({u'type': u'Point', u'coordinates': [-117.4365, 41.0826, 11.5]})
  • what exactly does the *v do in.. m.append(arcpy.PointGeometry(arcpy.Point(*v[2]), sr)) – ziggy Sep 29 '16 at 19:17
  • I know its some sort of unpacking generator but I always found it difficult to understand and incorperate – ziggy Sep 29 '16 at 19:17
  • 2
    You are correct, it is unpacking the points. So in the earthquakes list, the 3rd item in each item is a nested list [x, y]. By using the arcpy.Point(*list) syntax, it unpacks the list into the args for Point which are ({X}, {Y}, {Z}). So in this case, it is only unpacking the x and y coordinates to the points. So in short, arcpy.Point(*[x, y]) is the same as arcpy.Point(x, y). – crmackey Sep 29 '16 at 19:36
  • that would be super useful for polygon coordinates. neat trick – ziggy Sep 29 '16 at 20:16
  • Exactly, it is really nice when you are working with existing arrays of coordinates. – crmackey Sep 29 '16 at 20:40

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.