topic Re: Converting Shapefile to new coordinates in Python Questions

Converting Shapefile to new coordinates

deleted-user-Pi5x1UjkS2PY — Thu, 30 Oct 2014 14:33:13 GMT

I have the code below that is taking an SQL Query and making it into a feature class / shapefile that will be used in an JavaScript application. So it needs to be projected using WGS_1984_Web_Mercator. I know that the latitude/longitude in the table is not in this coordinate system.

When I import the shapefile created below into my project it is not projecting in WGS_1984_Web_Mercator as the points are not in the right spot. How can I fix this?

    cnxn = pyodbc.connect("DSN=AR")
    cursor = cnxn.cursor()
    cursor.execute("Select STORENAME, LATITUDE, LONGITUDE FROM GISWEB.VW_MCYWEB WHERE DISTRICTID = 949")

    datarray = []
    for row in cursor:
        datarray.append(tuple(row))

dtype = np.dtype([('STORENAME', '|S25'), ('xCoord', '<f8'), ('yCoord', '<f8')])
npyarr = np.array(datarray, dtype)

#webmercator 3857 - WGS_1984_Web_Mercator 3785
##SR = arcpy.Describe("C:/data/texas.gdb/fd").spatialReference

    sr = arcpy.SpatialReference(3785)
    #out_fc = "C:/inetpub/mystores3"
    out_fc = "//SP000XSSQL51/wwwroot/GIS maps and shapefiles/GISWeb/BaseLayers/mystorestest"
    if arcpy.Exists(out_fc):
        arcpy.Delete_management(out_fc)

arcpy.da.NumPyArrayToFeatureClass(npyarr, out_fc, ['xCoord','yCoord'], sr)

Re: Converting Shapefile to new coordinates

DanPatterson_Retired — Thu, 30 Oct 2014 14:42:22 GMT

setting it to the projection you want won't work you actually have to project the data...see the demonstration scripts for the Project Tool.

Re: Converting Shapefile to new coordinates

MelitaKennedy — Thu, 30 Oct 2014 15:59:10 GMT

What Dan said!

Also, you should use 3857, rather than 3785. It's a very messy story, but 3857 *is* the currently supported one and 3785 / 102113 has been deprecated (made obsolete).

Melita

Re: Converting Shapefile to new coordinates

deleted-user-Pi5x1UjkS2PY — Thu, 30 Oct 2014 18:44:31 GMT

I don't get what I am doing wrong here. I feel like the code I am using is not taking the actual points and putting them where they belong with the new coordinate system. I used the exact parameters another shapefile that I have that works in my MXD and used those to create a new shape file. The new shapefile shows those exact parameters in its projection file but all of points in the wrong place. You know ocean on top of each other unless you zoom in super tight. Any suggestions? Key code has been bolded.

def defAddSQLTable_CreateShapefile():
    cnxn = pyodbc.connect("DSN=AR")
    cursor = cnxn.cursor()
    cursor.execute("Select STORENAME, LATITUDE, LONGITUDE FROM GISWEB.VW_MCYWEB WHERE DISTRICTID = 949")

    datarray = []
    for row in cursor:
        datarray.append(tuple(row))

dtype = np.dtype([('STORENAME', '|S25'), ('xCoord', '<f8'), ('yCoord', '<f8')])
npyarr = np.array(datarray, dtype)

#webmercator 3857 - WGS_1984_Web_Mercator 3785
    sr = arcpy.SpatialReference(3857)
    #out_fc = "C:/inetpub/mystores3"
    out_fc = "//SP000XSSQL51/wwwroot/GIS maps and shapefiles/GISWeb/BaseLayers/mystorestest"
    if arcpy.Exists(out_fc):
        arcpy.Delete_management(out_fc)

arcpy.da.NumPyArrayToFeatureClass(npyarr, out_fc, ['xCoord','yCoord'], sr)

#Convert to Web Mercator

   # input data is in NAD 1983 UTM Zone 11N coordinate system
    input_features = r"//SP000XSSQL51/wwwroot/GIS maps and shapefiles/GISWeb/BaseLayers/mystorestest.shp"
    # output data
    output_feature_class = r"//SP000XSSQL51/wwwroot/GIS maps and shapefiles/GISWeb/BaseLayers/mystorestest2"
    # create a spatial reference object for the output coordinate system
    #out_coordinate_system = arcpy.SpatialReference(3857)
    out_coordinate_system = arcpy.SpatialReference(102100)
    # run the tool

#arcpy.Project_management(input_features, output_feature_class, out_coordinate_system)

### Copied form another layer that is working using my old non python way.
    arcpy.Project_management(input_features, output_feature_class, \
                                     r"PROJCS['WGS_1984_Web_Mercator'," + \
                                     r"GEOGCS['GCS_WGS_1984_Major_Auxiliary_Sphere'," + \
                                     r"DATUM['D_WGS_1984_Major_Auxiliary_Sphere'," + \
                                     r"SPHEROID['WGS_1984_Major_Auxiliary_Sphere',6378137.0,0.0]]," + \
                                     r"PRIMEM['Greenwich',0.0]," + \
                                     r"UNIT['Degree',0.0174532925199433]]," + \
                                     r"PROJECTION['Mercator']," + \
                                     r"PARAMETER['False_Easting',0.0]," + \
                                     r"PARAMETER['False_Northing',0.0]," + \
                                     r"PARAMETER['Central_Meridian',0.0]," + \
                                     r"PARAMETER['Standard_Parallel_1',0.0]," + \
                                     r"UNIT['Meter',1.0]]")

Re: Converting Shapefile to new coordinates

DanPatterson_Retired — Thu, 30 Oct 2014 18:49:12 GMT

without examining the code...when you created the points, ensure that you specified longitude then latitude (ie X first, then Y) in that order...I have had two coordinate switcheroo cases this morning already.

Re: Converting Shapefile to new coordinates

MelitaKennedy — Fri, 31 Oct 2014 18:02:04 GMT

Are you positive that this shapefile is also using UTM 11N? If that's the defined coordinate system, try to verify that it's correct. Add it to ArcMap and change the display units to decimal degrees--do the values make sense?

Melita