1

I'm creating a raster grid in ArcPy and would like to iterate over it and get each cell's center. The raster grid is in UTM but I'm issuing the coordinates in as GCS WGS 1984 XY. Here's what I'm doing:

import arcpy
import os


wgs84 = arcpy.SpatialReference(4326)
utm = arcpy.SpatialReference(32617)

lower_left = arcpy.PointGeometry(arcpy.Point(X=82.4819495, Y=34.2706214), spatial_reference=wgs84).projectAs(utm)
upper_right = arcpy.PointGeometry(arcpy.Point(X=82.0848611, Y=34.5154197), spatial_reference=wgs84).projectAs(utm)
grid_extent = arcpy.Extent(
    lower_left[0].X,
    lower_left[0].Y,
    upper_right[0].X,
    upper_right[0].Y,
    0, 0, 0, 0,
    utm
)

cell_size = 25
pixeltype = "F32"

# make the rasterinfo object
ras_info = arcpy.RasterInfo()
ras_info.setSpatialExtent(utm)
ras_info.setExtent(grid_extent)
ras_info.setCellSize((cell_size, cell_size))
ras_info.setPixelType(pixeltype)

# make the raster object
ras = arcpy.Raster(ras_info)

# Now convert it to points to get all the cell centers
ras_cell_centers = arcpy.RasterToPoint_conversion(ras, os.path.join(r"in_memory", "rcc"))

But, when I do

with arcpy.da.SearchCursor(ras_cell_centers, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        print(row[0])

it gives me no point objects (but no errors either).

Where am I going wrong? I would like any solutions that are not IO bound and allows me to do this in memory. I checked out the Raster Cell Iterator and can use it, but it only gives me column-row not the cell's center projected coordinates.

1 Answer 1

1

There are multiple issues with the code:

  1. You appeared to be mixing ArcMap and ArcPro flavours of arcpy so this code would probably have failed.
  2. There is no property called setSpatialExtent() on ras_info so this code would never have worked.
  3. The ProjectAs() on your point geometry was returning a None object, seems like a bug in arcpy, so you guessed it this code would never have worked....
  4. Your logic is flawed, you "create" a raster but never assign any cell values so the conversion to points is turning NODATA into nothing, so it would never have worked! :)

To overcome the issue with ProjectAs() failing I determined the UTM equivalent coordinate on epsg.io.

So I was able to get your code working as shown below:

import arcpy
arcpy.env.addOutputsToMap = True

utm = arcpy.SpatialReference(32617)
ll = arcpy.Point(X=2028185.68, Y=16077634.07) # X=82.4819495, Y=34.2706214
ur = arcpy.Point(X=2060559.50 , Y=16043794.27 ) # X=82.0848611, Y=34.5154197
lower_left = arcpy.PointGeometry(ll, spatial_reference=utm)
upper_right = arcpy.PointGeometry(ur, spatial_reference=utm)

# create extent
grid_extent = arcpy.Extent(lower_left.centroid.X, lower_left.centroid.Y, upper_right.centroid.X, upper_right.centroid.Y, 0, 0, 0, 0,utm)

cell_size = 25
pixeltype = "F32"

# make the rasterinfo object
ras_info = arcpy.RasterInfo()
ras_info.setSpatialReference(utm)
ras_info.setExtent(grid_extent)
ras_info.setCellSize((cell_size, cell_size))
ras_info.setPixelType(pixeltype)

# make the raster object and assign cell values 1
ras = arcpy.Raster(ras_info)
for i,j in ras:
    ras[i,j] = 1

# Now convert it to points to get all the cell centers
ras_cell_centers = arcpy.RasterToPoint_conversion(ras, r"memory\rcc")

# print coordinates
with arcpy.da.SearchCursor(ras_cell_centers, ["SHAPE@XY"]) as cursor:
    for row in cursor:
        print(row[0])

WARNING: your extent and cell size generates a raster of ~ 1.7 million cells so the print XY coordinates takes a long time to complete

Finally here is you logic crushed down to 5 lines of code and I suspect runs considerably faster.

import arcpy
arcpy.env.outputCoordinateSystem = arcpy.SpatialReference(32617)
out_raster = arcpy.sa.CreateConstantRaster(1, "INTEGER", 25, "2028198.18 16043806.77 2060548.18 16077631.77")
out_raster.save(r"memory\test")
ras_cell_centers = arcpy.RasterToPoint_conversion(out_raster, r"memory\rcc2")
1
  • 1
    Thanks! For #1: Was using the Pro documentation. For #2: That was a mis-copy on my part as I am writing this on an intranet where I can’t do a copy-paste. For #3: For the end result, users will be issuing coordinates in GCS lat lon so I’ll work on getting that functional. For #4 you’re right, I am new to dealing with raster data so your answer was very helpful to me.
    – auslander
    Apr 27 at 22:27

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.