1 Reply Latest reply on May 11, 2012 8:29 AM by davisam1

    calculate nearest ag pixel for each urban pixel in a raster

    davisam1
      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
        • Re: calculate nearest ag pixel for each urban pixel in a raster
          davisam1
          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 = "XX_" + 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, "")