+ Reply to Thread
Results 1 to 2 of 2

Thread: calculate nearest ag pixel for each urban pixel in a raster

  1. #1
    Amelie Davis
    Join Date
    Aug 2011
    Posts
    22
    Points
    0
    Answers Provided
    0


    0

    Default calculate nearest ag pixel for each urban pixel in a raster

    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

  2. #2
    Amelie Davis
    Join Date
    Aug 2011
    Posts
    22
    Points
    0
    Answers Provided
    0


    0

    Default Re: calculate nearest ag pixel for each urban pixel in a raster

    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.

    Code:
    # ---------------------------------------------------------------------------
    # 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, "")

+ Reply to Thread

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts