1

I have a set of 100 points and 100 rasters. Points are stored in one shapefile over large spatial extent and small, not-overlapping rasters are stored in gdb. Rasters are something like tiles. Together rasters cover the whole extent of points, but one by one, only one raster by one point can overlap. That is, many points are outside of particular rasters.

Here is an example: each point has its own one or two rasters, but rasters are not stored in mosaic, rather one by one:

Example: each point has its own one or more rasters, but rasters are not continuous as mosaic

How do I Extract raster (multiple) values to points for each raster to overlapping point?

Normally, I would firstly merge the individual rasters into new mosaic (Mosaic to new raster) and then simply Extract (Multi) Values to Points. NOw, I want to avoid the mosaic step due to memory issues (mosaic creates raster of 36 GB).

I was thinking of doing something like this, using Python and ArcPy:

Pseudocode:

1. Read points shapefile
2. Read rasters into raster list
3. Loop through each row of shapefile (`arcpy.da.seach.coursor`)
4. Check if the shapefile row overlays one of rasters in raster lists
    4.1. if yes - extract raster value to point, move to next row
    4.2. in no - move to next row

However, I am not sure how to do this overlay check.

Here is where I am now:

import arcpy, os
arcpy.CheckOutExtension("spatial")

# Define path to working directory
inWD = "C:/Users/extract" 

# set working environment
arcpy.env.workspace = os.path.join(inWD, "analyzed.gdb")
arcpy.env.overwriteOutput = True

# Define input data
inPoints =  os.path.join(inWD,  "pts.shp") 

# Read all raster in a list
rasters = arcpy.ListRasters("rst*", "GRID")

# Loop throught points row by row
with arcpy.da.SearchCursor(inPoints ,'OID@') as cursor:
    for row in cursor:
        # Create a query to subset rows one by one
        sql="""{0} = {1}""".format(arcpy.AddFieldDelimiters(inPoints ,
                                                            arcpy.Describe(inPoints ).OIDFieldName),row[0])
        print sql

        # Create new feature layer row by row
        outLyrName = "lyr_" + str(row[0])
        arcpy.MakeFeatureLayer_management(inPoints, outLyrName, sql)

        # Check if point overlays one of rasters:
        !! HOw to check this??

        # Process: 
        arcpy.sa.ExtractMultiValuesToPoints(outLyrName, rasters, )

Simply extracting all points from all rasters in a list creates columns of empty data, so I think I can somehow do it one by one.

0

1 Answer 1

4

No need to loop through each row in a points shapefile.

Create a Mosaic dataset and add all rasters to it. The Mosaic dataset will read the rasters where they are stored on disk so the Mosaic dataset itself will be very small. Using catalog create a file geodatabase then right click it and:

enter image description here

The output Mosaic dataset can be then easily used within Extract Multi Values to Points, without creating the empty fields.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.