4

I am trying to extract points from multiple raster files into a points shapefile. I have about 320,000 unique XY coordinates which need to extracted values from 37 raster files. I tried running the following code and it has been three days since it has started.

import arcpy
import os
from arcpy.sa import *

if arcpy.CheckExtension("spatial") == "Unavailable":
    sys.exit("No ArcGIS Spatial Analyst licence available - exiting")
else:
     arcpy.CheckOutExtension("spatial")
#     arcpy.CheckInExtension("spatial")

arcpy.env.workspace = r"H:\GIS Project\in_rasters"
arcpy.env.overwriteOutput = True
arcpy.env.parallelProcessingFactor = "50"
point_feature_class = r"H:\GIS Project\XY_Points.shp"
rasters = arcpy.ListRasters("*","TIF")


for raster in rasters:
    ExtractMultiValuesToPoints(point_feature_class, raster, 'NONE')

My computer has enough ram (16 gigs) and I am running it in parallel. I would think this function would take a couple of hours maximum based off this post Getting raster values for large number of point features?. Albeit, this post is not exactly the same, it still shows some anecdotal evidence for processing time with large data.

My raster files are all projected the same (WGS 84) however they do not have the same resolution.

Is it normal for this function to take a long time to finish?

  • How long does one raster (typical of your rasters) take to process? Have you tried the solution posted at gis.stackexchange.com/a/94301/115? – PolyGeo Jul 3 '18 at 0:06
  • That is my next step is to try one raster layer at a time. Sadly I am sharing a workstation and wont have access for a few hours. The solution posted at gis.stackexchange.com/a/94301/115 is similar to the approach I have here. – Douglas Jul 3 '18 at 0:12
  • 1
    This is not an area in which I specialize but I wonder whether the parallel processing factor is causing the processes to "trip over each other". I wonder whether running them sequentially may be better in this instance. – PolyGeo Jul 3 '18 at 0:15
  • Have you tried adding a spatial index to your shapefile? This will speed things up, it certainly did for me when testing 320,000 random points against a raster. It took about 20 minutes to sample 3 rasters (all in one go) with 320,000 points. Not hours! – Hornbydd Jul 3 '18 at 10:19
2

According to the help on this tool the environment parallel processing factor isn't considered.. not many tools do parallel process.

You might get a performance increase by running all rasters at the same time but on a smaller chunk of points:

import arcpy
import os
from arcpy.sa import *

if arcpy.CheckExtension("spatial") == "Unavailable":
    sys.exit("No ArcGIS Spatial Analyst licence available - exiting")
else:
     arcpy.CheckOutExtension("spatial")
#     arcpy.CheckInExtension("spatial")

arcpy.env.workspace = r"H:\GIS Project\in_rasters"
arcpy.env.overwriteOutput = True
arcpy.env.parallelProcessingFactor = "50" # not used by this tool
point_feature_class = r"H:\GIS Project\XY_Points.shp"
rasters = arcpy.ListRasters("*","TIF")

# build a list of your rasters
AllRasters = [] # an empty list
for raster in rasters:
    fName, fExt = os.path.splitext(raster)
    AllRasters.append([raster,fName])

FeatCount = int(arcpy.GetCount_management(point_feature_class).getOutput(0))
MemFCList = []
ChunkSize = 1000 # how many features to do at the same time
CountList = range(0,FeatCount,ChunkSize)

for StartValue in CountList:
    arcpy.AddMessage('Running chunk {}'.format(StartValue))
    MemFC = 'in_memory\pfc_{}'.format(StartValue)
    MemFCList.append(MemFC)
    # select a chunk of data
    arcpy.Select_analysis(point_feature_class,MemFC,'FID >= {} AND FID < {}'.format(StartValue,StartValue + ChunkSize))
    # run the tool on all rasters at once.
    ExtractMultiValuesToPoints(MemFC, AllRasters, 'NONE')

# replace the original by merging the chunks
arcpy.Merge_management(MemFCList,point_feature_class)

Each chunk of chunksize (default 1000) is copied into the in_memory workspace which should speed things up if you're accessing the features from a slow workspace (either a slow drive or network storage) though it will not help if your rasters are on a slow drive or network storage.

A few things to check first:

  1. Repair geometry on your input points, there could be some dud points that are gumming up the works.
  2. Ensure your rasters are on a local, preferable fast, drive. Accessing rasters on slow laptop drives, USB 2 and network/cloud drives will all make this process tediously slow.
  3. The rasters are uncompressed. Having to decompress your rasters constantly uses CPU cycles that are better used elsewhere.
  4. Your rasters aren't in a slow format like ASCII or XYZ. Both of these are rubbish for processing, if you do have either of these formats then converting to GeoTIFF or ERDAS IMG will reduce your processing time notably.
  • all of my rasters are in TIFF format, so I think I am okay. I am working on a network drive, so I will move over all my files to a local drive to see if it fixes the speed problem. I am currently waiting for a spatial analyst package to be checked back in. I will let you know my results later! – Douglas Jul 3 '18 at 4:26
  • I think dividing the area, say via fishnet, and using the extent of the each cell as processing extent might be a better approach rather than using FID ranging (which is a bit unreliable without knowing the FID field range). – fatih_dur Jul 3 '18 at 5:10
  • @fatih_dur the range of FID is 0 to count of features minus 1 except when you're editing the shapefile. It's one of the handy quirks in a shapefile that almost make up for the common shortfalls. Caveat: this does not apply to any other feature storage type (file/personal geodatabase for example), these feature storage types get an OBJECTID assigned at creation and keep the same ID for life. – Michael Stimson Jul 3 '18 at 5:32
  • @MichaelStimson, knowledge of the day to me then, cheers. – fatih_dur Jul 3 '18 at 5:33
  • 1
    @MichaelStimson so something like this at the end of the code? unique_name = arcpy.CreateUniqueName("temp.shp") arcpy.Merge_management(MemFCList,unique_name) arcypy.Rename_management (unique_name, locust_data.shp) – Douglas Jul 3 '18 at 22:58

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.