calculate nearest ag pixel for each urban pixel in a raster

Question

Hi!
I have a raster (NLCD2006) and I would like to calculate the distance of the nearest pixel of one class from another class.  
I am struggling to think through what would be the best approach.  
I thought about reclassifying all the ag and urban pixels to one, setting the rest to NoData and running the euclidean distance tool but that isn't quite what I want since it might give me the nearest a pixel rather than the urban one.  
Then I thought I could turn the ag pixels to points, same for the urban pixels but in a separate shapefile and then run the near tool but I have too many cells/points for this approach (I'm running this for the entire US for 30 by 30m cells.  
I can't think of another approach?  
Does anyone have any suggestions?  
Any help is greatly appreciated,
Amelie

AmelieDavis · Answer

OK I think I figured out how to do it so I figure I would post the answer (even if I did start this thread...)
I'm not an expert coder so there probably is a much more elegant/faster way to do this but at least it works.  
Hopefully this helps someone.  
# ---------------------------------------------------------------------------
# Dist2LandCover_v1.py
# Created on: 2012-05-09
# By Amelie Davis 
# Description: 
# ---------------------------------------------------------------------------

############ IMPORTANT ########################
# Assumes that your land cover layer is projected AND cell size is 30 meters
# Uses NLCD 2006 classification codes
# output file must be a shapefile and the output name must have .shp at the end of it.
##############################################

# Set the necessary product code
# import arcinfo

# Import arcpy module
import arcpy
from arcpy import env

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

###################### Script arguments

Input_Fishnet_with_unique_ID2_field = arcpy.GetParameterAsText(0)
#Input_Fishnet_with_unique_ID2_field = "D:\DATA\USA\MigBirdFishnet\NLCDfishnet20mi_p.shp" # provide a default value if unspecified

Input_Land_Cover__NLCD_2006_ = arcpy.GetParameterAsText(1)
#Input_Land_Cover__NLCD_2006_ = "D:\DATA\NLCD\nlcd2006_cook" # provide a default value if unspecified

Output_Workspace = arcpy.GetParameterAsText(2)
#Output_Workspace = "D:\DATA\DELETE" # provide a default value if unspecified

# Set environment settings
env.workspace = Output_Workspace
env.mask = Input_Land_Cover__NLCD_2006_

Zone_value_you_want_to_calc_dist_from = arcpy.GetParameterAsText(3)
#Zone_value_you_want_to_calc_dist_from = "11" # provide a default value if unspecified

LC_codes = arcpy.GetParameterAsText(4)

Final_Output_shp = Output_Workspace + "\" + arcpy.GetParameterAsText(5)

####################

# Process: Reclassify #1
#set to 1 LC code from which you want to calculate distance.  Example: water = "11"
text_in = Zone_value_you_want_to_calc_dist_from+" "+Zone_value_you_want_to_calc_dist_from+" 1"
fromClassremapList = text_in.split()
fromClassremapList = [int(k) for k in fromClassremapList]
remap = arcpy.sa.RemapRange([fromClassremapList]) 
nlcd_wat = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap, "NODATA")
#remap = arcpy.sa.RemapRange([[11,11,1],[12,95,"NODATA"]]) 
#nlcd_wat = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap, "DATA")
out_rc1 = Output_Workspace + "\out1"
nlcd_wat.save(out_rc1)

# Process: Euclidean Distance
out_euc = arcpy.sa.EucDistance(nlcd_wat, "", "30", "")
out_euc_s = Output_Workspace + "\euc_dist_1"
out_euc.save(out_euc_s)
arcpy.AddMessage("finished calculating euclidian distance")

# Process: Int
out_Int = arcpy.sa.Int(out_euc)
out_Int_s = Output_Workspace + "\water_dist"
out_Int.save(out_Int_s)

# Process: Copy initial fishnet shapefile so that original is not modified
arcpy.Copy_management(Input_Fishnet_with_unique_ID2_field, Final_Output_shp, "")

# Process: Reclassify (2)
#set to 1 LC code for which you want to know the average distance from LC in Process: Reclassify #1
#remap2 = arcpy.sa.RemapRange([[11,94,0],[95,95,1]]) #old
#remap2 = arcpy.sa.RemapRange([[95,95,1]]) #works
LC_CODE = LC_codes.split()
for i in LC_CODE:
    text_in2 = i+" "+i+" 1"
    forClassremapList = text_in2.split()
    forClassremapList = [int(k) for k in forClassremapList]
    remap2 = arcpy.sa.RemapRange([forClassremapList]) 
    nlcd_ag = arcpy.sa.Reclassify(Input_Land_Cover__NLCD_2006_, "VALUE", remap2, "NODATA")
    out_rc2 = Output_Workspace + "\out2"
    nlcd_ag.save(out_rc2)
    arcpy.AddMessage("finished reclassifying")

# Process: Times
    # save output permanently?
    out_times = arcpy.sa.Times(out_Int, nlcd_ag)
    out_temp2 = Output_Workspace + "\temp_ag"
    out_times.save(out_temp2)

field_name = "xyz_" + str(i)

# Process: Add field to initial fishnet shapefile
    arcpy.AddField_management(Final_Output_shp, field_name, "SHORT", 6, "", "","refcode", "NULLABLE")

# Process: Calculate Field
    arcpy.CalculateField_management(Final_Output_shp, field_name, "!FID!","PYTHON_9.3","")

arcpy.AddMessage("Starting zonal statistics")
    # Process: Zonal Statistics as Table
    outTable = Output_Workspace + "\dist_per_cell.dbf"
    outZSaT = arcpy.sa.ZonalStatisticsAsTable(Final_Output_shp, field_name, out_times, outTable, "DATA", "ALL")
    arcpy.AddMessage("finished calculating mean +/- std distance per cell")

# Process: Join Field
    fieldList = ["MEAN", "STD", "MIN"]
    arcpy.JoinField_management(Final_Output_shp, field_name, outTable, field_name,fieldList)

arcpy.AddMessage("finished with land cover" + str(i))

# Delete intermediary files
##arcpy.Delete_management(out_rc1, "")
##arcpy.Delete_management(out_Int_s, "")
##arcpy.Delete_management(out_rc2, "")
##arcpy.Delete_management(out_euc_s, "")