Latitude & Longitude rasters in ArcGIS 10

The ability to calculate latitude and longitude rasters using the GRID commands $$YMap and $$XMap was removed in ArcGIS 10 (as well as all the GRID commands including $$COLMAP, $$ROWMAP, $$NROWS, $$NCOLS, $$CELLSIZE). To get around this you can still load the ArcGIS 9.3 geoprocessor in ArcGIS 10 and then call the functions using the Single Output Map Algebra tool. Latitude & longitude rasters can also be calculated with NumPy or the FlowAccumulation function but these aren't great options for really large rasters. Make sure to use a DEM that is in a geographic coordinates system and then project the results.

def lat_long_rasters(ned_raster, lat_raster, lon_raster):
    import arcgisscripting
    gp = arcgisscripting.create(9.3)
    gp.CheckOutExtension("Spatial")
    gp.WorkSpace = os.path.dirname(ned_raster)
    gp.SingleOutputMapAlgebra_sa(
        "Con(IsNull({0}), {1}, $$YMap)".format(ned_raster, ned_raster), 
        lat_raster)
    gp.SingleOutputMapAlgebra_sa(
        "Con(IsNull({0}), {1}, $$XMap)".format(ned_raster, ned_raster), 
        lon_raster)
    del gp