0

I am trying to update the area field in sqkm of my table using Search and Update cursor from the arcpy library. The area needs to be updated using the different utm zone. I intersected my shp (that I called "outIntersectPath") with the utm zone shp. Now my intersected shp has a column with the UTM ZONE that goes from 1 to 60 (integer) and another column NTH_STH that represents the Northern ("N") and the Southern ("S") hemisphere.

This is the code:

def utmConversion(outIntersectPath):
 # calculate Area per UTM Zones
 count = 1
 utmNorthZone = 32600
 utmSouthZone = 32700
 field = ["SHAPE@AREA", "BASIN_AREA_SQKM", "UTM_ZONES", "NTH_STH"]
 arcpy.AddField_management(outIntersectPath, "BASIN_AREA_SQKM", "DOUBLE")
 print "Field {0} has been added to the {1}".format(field[1], outIntersectPath)

 while count < 61:

    with arcpy.da.SearchCursor(outIntersectPath, field) as cursor:
        for row in cursor:
            utmZone = row[2]
            hemisphere = row[3]
            if utmZone == count and hemisphere == "N":
                utmNorthZone += 1
                count += 1
                # print utmZone
                # print utmNorthZone
                out_coordinate_system = arcpy.SpatialReference(utmNorthZone)
                print out_coordinate_system.name

                with arcpy.da.UpdateCursor(outIntersectPath, field, "", out_coordinate_system) as ncursor:
                    for nrow in ncursor:
                        nrow[1] = nrow[0] * 10 ** 6
                        ncursor.updateRow(nrow)
                print "Field {0} has been updated to {1}.".format(field[1], outIntersectPath)
            elif utmZone == count and hemisphere == "S":
                utmSouthZone += 1
                count += 1
                out_coordinate_system = arcpy.SpatialReference(utmSouthZone)
                print out_coordinate_system.name
                with arcpy.da.UpdateCursor(outIntersectPath, field, "", out_coordinate_system) as scursor:
                    for srow in scursor:
                        srow[1] = srow[0] * 10 ** 6
                        scursor.updateRow(srow)
                print "Field {0} has been added to {1}.".format(field[1], outIntersectPath)
            else:
                print "The Process failed..."

What I am doing wrong?.. it is not calculating the area per each UTM Zone.

  • why is your calculation nrow[0] * 10 ** 6, the geometry area should be in square metres which means divide by 10**6. What spatial reference is the data in natively? I think you are misunderstanding the spatial reference parameter of the cursor but I'm fairly sure I can find a workaround. – Michael Stimson Jul 20 '17 at 21:05
  • Yeah I wanted to divide to get the sq km. That is incorrect.. the native proj is wgs84. – ManuCiao Jul 20 '17 at 21:41
  • WGS84 is compatible with UTM as specified in your script. You can run into problems in projection if the source spatial reference has an incompatible datum to the projected spatial reference, if you encounter this situation (like NAD->UTM I think, certainly AGD66->WGS84 UTM) it would require a transformation parameter in the projectAs() line. – Michael Stimson Jul 20 '17 at 21:44
2

Your script is fairly well written but is unnecessarily complicated and I think you have the wrong idea of what the spatial_reference parameter of the update cursor is for, from the online help:

The Spatial Reference of the feature class can be specified with either a SpatialReference object or string equivalent.

There is no indication that geometries from the cursor will be projected.

I can't really follow your loop logic, using a while, search loop and update loop; in the end all features would be updated to the zone in question as you're omitting the where_clause in the update cursor. I think it would be best to just loop on one update cursor and decide whether or not to update the feature and what zone rather than breaking down into branches with almost duplicated code where all that needs to change is the spatial reference.

I think this is what you're after:

def utmConversion(outIntersectPath):
    utmNorthZone = 32600
    utmSouthZone = 32700
    # note I've changed SHAPE@AREA to SHAPE@ because I want a Geometry object
    # to project and then return the area of rather than a double.
    field = ["SHAPE@", "BASIN_AREA_SQKM", "UTM_ZONES", "NTH_STH"] 

    if arcpy.ListFields(outIntersectPath,"BASIN_AREA_SQKM") == None:
        # add the intended field if it doesn't exist
        print "Adding Field BASIN_AREA_SQKM to the {}".format(outIntersectPath)
        arcpy.AddField_management(outIntersectPath, "BASIN_AREA_SQKM", "DOUBLE")

    with arcpy.da.UpdateCursor(outIntersectPath, field) as UCur:
        for URow in UCur:
            Geom    = URow[0]
            UTMZone = URow[2]
            Hemi    = URow[3].upper() # always compare characters in either upper or lower.
            UTMZone += 1 # not sure why you're doing this but I'll assume it's for a good reason.

            if UTMZone > 0 and UTMZone < 61:
                # only calculate and update where the zone is in the expected range
                if Hemi == 'N':
                    OutSR = arcpy.SpatialReference(utmNorthZone + UTMZone)
                elif Hemi == 'S':
                    OutSR = arcpy.SpatialReference(utmSouthZone + UTMZone)
                else:
                    OutSR = None # not properly defined hemisphere

                if OutSR != None:
                    # Spatial reference is properly defined..
                    # only calculate and update rows where hemisphere is
                    # properly defined
                    ProjectedGeom = Geom.projectAs(OutSR)
                    AreaSqKm = ProjectedGeom.area / (10**6) # turn square metres into square kilometres.
                    URow[1]  = AreaSqKm
                    UCur.updateRow(URow)

I can't see why you're incrementing the UTM zone by one but can only assume that the values in the field are reduced by one and you're compensating for it.

Like I said in my comment, the area will be returned in square metres, to convert to square kilometres you must divide not multiply, so the calc is AreaSqKm = Geom.area / (10**6) to achieve this I grab the geometry object of the feature, project and then get the area property. This way you're projecting on-the-fly for each feature as encountered, provided the zone and hemisphere are populated correctly.

  • 1
    The only thing I need to change it is: proj = Geom.projectAs(OutSR) AreaSqKm = proj.area / (10**6) – ManuCiao Jul 21 '17 at 12:47
  • ..otherwise it won't calculate the area of the projected geom obj. – ManuCiao Jul 21 '17 at 13:51
  • @ManuelaSabatino Ah yes, you're right, thank you for pointing that out. The python projection (projectAs) returns a geometry in the new spatial reference... I must have had a flashback to ArcObjects where the Geom.project(SpatialReference) has no return value but replaces the geometry object with a new object in the desired spatial reference. I'll change that right now. – Michael Stimson Jul 23 '17 at 21:22

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.