Memotong raster dengan beberapa dataset atau poligon?

8

Saya ingin klip DEM menggunakan kotak poligon. Mungkin lebih mudah untuk menggunakan banyak poligon dalam satu file bentuk, tetapi saya belum berhasil ini jadi saya mencoba menggunakan for for loop sehingga saya bisa mengulang setiap dataset dalam gdb (masing-masing hanya berisi satu poligon).

Ini kode saya (melakukannya di jendela python).

#creating a workspace and a list of feature classes
arcpy.env.workspace = "C:/data/lidar/lidar.gdb"
fcs = arcpy.ListFeatureClasses()

#looping through each feature class and creating a raster based on the extent of   
#feature class

for fc in fcs:
    arcpy.Clip_management("perth", "#", "C:/data/lidar", fc, "", "ClippingGeometry")

Namun kode saya tidak dieksekusi, hanya duduk di sana, menunggu sesuatu yang lain ... tapi apa? Saya bisa membuatnya bekerja untuk satu klip, tetapi tidak dengan loop.

Saya yakin saya harus melakukan sesuatu yang lain untuk output, untuk memberi nama setiap raster baru berdasarkan kelas fitur atau sesuatu ... tapi sekali lagi, tidak tahu caranya. Harap beri tahu saya jika saya harus menambahkan info lagi.

Rosie Bell
sumber
2
Mungkin uji terlebih dahulu apakah loop Anda berfungsi dengan mengomentari Klip dan menempatkan cetakan atau AddMessage untuk melihat apakah itu memuntahkan nama setiap kelas fitur OK.
PolyGeo
Bisakah Anda menjelaskan mengapa Anda perlu melakukan kliping ini? Kecuali jika output diperlukan untuk komunikasi dengan perangkat lunak lain, biasanya ada metode yang jauh lebih efisien untuk menganalisis data raster poligon-demi-poligon, jadi mungkin jawaban terbaik untuk pertanyaan Anda adalah dengan menyarankan pendekatan yang berbeda sama sekali.
whuber
Permintaan maaf untuk memposting dan menjalankan - terima kasih semuanya atas bantuan Anda dan saya akan kembali untuk menindaklanjuti ini, sesegera mungkin.
Rosie Bell
wuber, saya ingin melakukan ini sehingga saya bisa mengiris DEM kami menjadi potongan yang bisa diatur untuk kemudian membuat kontur dan kemudian menjahit kembali bersama. Mungkin bukan cara paling efisien untuk melakukannya.
Rosie Bell
1
Perangkat lunak apa yang digunakan? ArcGIS, QGIS, dll?
Chad Cooper

Jawaban:

6

Satu hal yang saya perhatikan adalah bahwa parameter ketiga Anda adalah output kode keras (C: / data / lidar). Cara ini ditulis sekarang akan mengulangi masing-masing fitur Anda dan menimpa output setiap kali, tetapi karena Anda mungkin tidak mengizinkan penimpaan file secara otomatis, ini bisa berpotensi menjadi hang-up. Coba buat nama keluaran unik untuk setiap iterasi:

#creating a workspace and a list of feature classes
arcpy.env.workspace = "C:/data/lidar/lidar.gdb"
fcs = arcpy.ListFeatureClasses()

#looping through each feature class and creating a raster based on the extent of   
#feature class

i=0
for fc in fcs:
    outputPath = "C:/data/lidar" + str(i)
    i+=1
    arcpy.Clip_management("perth", "#", outputPath, fc, "", "ClippingGeometry")

Juga, saya tidak yakin bahwa Anda bermaksud untuk menempatkan output di folder C: / data bernama lidar ... perhatikan bahwa parameter ketiga dalam klip adalah path lengkap dari raster output Anda, bukan folder yang akan ditempatkan in. Jika Anda tidak menentukan ekstensi dalam nama jalur keluaran Anda dan menempatkannya di folder standar, itu akan berupa kisi, jadi saat ini program Anda sedang mencoba membuat dataset kisi baru bernama 'lidar' di C: / folder data.

kaki biru
sumber
Banyak terima kasih mbenedetti, itu adalah kurangnya nama yang unik untuk setiap klip baru yang merupakan salah satu masalah saya. Saya juga mengalami sesuatu yang aneh dengan itu menjadi gdb - ketika saya mencoba menggunakan shapefile di folder, itu berfungsi dengan baik. Jadi saya mungkin tidak memiliki jalur file yang benar (tidak terbiasa bekerja dengan gdb).
Rosie Bell
5

untuk pencari masa depan: Berikut adalah versi modifikasi dari skrip alat perpecahan USGS raster yang tidak memerlukan apa pun di atas tingkat lisensi ArcGIS Basic (ArcView):

"""
Raster Split Tool 6/16/2011
ArcGIS 10 Script Tool
Python 2.6.5

Contact:
Douglas A. Olsen
Geographer
Upper Midwest Environmental Sciences Center
U.S. Geological Survey
2630 Fanta Reed Road
La Crosse, Wisconsin 54603
Phone: 608.781.6333

#####
modified 2015-04-15 by Jeremiah Poling ([email protected])
Now using only tools available at the ArcGIS Basic license level 
#####

"""
import arcpy
from arcpy import env
import os

# Get Split Shapfile Name
splitShape = arcpy.GetParameterAsText(0)
# Get Field Name
splitField = arcpy.GetParameterAsText(1)
# Get Output Directory
outDirectory = arcpy.GetParameterAsText(2)
# Get Raster to Split
splitRaster = arcpy.GetParameterAsText(3)
# Get type of output tif, img, or GRID
rasterType = arcpy.GetParameterAsText(4)
# Set Workspace environment
env.workspace = outDirectory

# Loop through the rows in the clipping shapefile
cursor = arcpy.SearchCursor(splitShape)
for row in cursor:
    resultName = row.getValue(splitField)

    # Create feature layer of current clipping polygon
    whereClause = '"' + splitField + '" = ' + "'" + resultName + "'"
    arcpy.MakeFeatureLayer_management(splitShape, 'currentMask', whereClause)

    # Replace spaces with underscores
    resultName = resultName.replace(' ','_')

    # Remove .shp suffix
    if resultName[-4:] == '.shp':
        resultName = resultName[:-4]

    arcpy.AddMessage("Processing: " + resultName)

    if rasterType == 'img':
        resultName = resultName + "." + rasterType
    elif rasterType == 'tif':
        resultName = resultName + "." + rasterType
    else:
        print ('No extension')

    # Save the clipped raster
    arcpy.Clip_management(
        in_raster = splitRaster,
        rectangle = "#",
        out_raster = resultName,
        in_template_dataset = 'currentMask',
        nodata_value="255",
        clipping_geometry="ClippingGeometry",
        maintain_clipping_extent="NO_MAINTAIN_EXTENT"
        )

    if arcpy.Exists('currentMask'):
        arcpy.Delete_management('currentMask')
Jeremiah Poling
sumber
4

Beberapa ide:

  1. Coba Raster Split Tool , tersedia sebagai alat skrip dari USGS (lihat kode sumber terlampir).
  2. Untuk kliping / ubin raster sederhana, gunakan alat ArcGIS 10.1 bawaan yang disebut Split Raster. Anda dapat membagi raster dengan jumlah ubin atau ukuran ubin.

masukkan deskripsi gambar di sini

Kode sumber untuk alat USGS Raster Split:

"""
Raster Split Tool 6/16/2011
ArcGIS 10 Script Tool
Python 2.6.5

Contact:
Douglas A. Olsen
Geographer
Upper Midwest Environmental Sciences Center
U.S. Geological Survey
2630 Fanta Reed Road
La Crosse, Wisconsin 54603
Phone: 608.781.6333
"""
import arcpy
from arcpy import env
import os
from arcpy.sa import *

# Get Shapefile Name
inShape = arcpy.GetParameterAsText(0)
# Get Split Shapfile Name
splitShape = arcpy.GetParameterAsText(1)
# Get Field Name
splitField = arcpy.GetParameterAsText(2)
# Get Output Directory
outDirectory = arcpy.GetParameterAsText(3)
# Get Raster to Split
splitRaster = arcpy.GetParameterAsText(4)
# Get type of output tif, img, or GRID
rasterType = arcpy.GetParameterAsText(5)
# Set Workspace environment
env.workspace = outDirectory

# Run Split command on Input Shapefile
arcpy.Split_analysis(inShape, splitShape, splitField, outDirectory)

# Loop through a list of feature classes in the workspace
for fc in arcpy.ListFeatureClasses():

    # Execute ExtractByMask
    rfc = ExtractByMask(splitRaster, fc)

    # Clean up shp file from work directory
    arcpy.Delete_management(fc, "")    

    arcpy.AddMessage("Processing: " + fc)

    # Replace spaces with underscores
    fc = fc.replace(' ','_')

    # Remove .shp suffix
    fc = fc[:-4]

    # Trim to 13 chars
    fc = fc[0:13]

    if rasterType == 'img':
        fc = fc + "." + rasterType
    elif rasterType == 'tif':
        fc = fc + "." + rasterType
    else:
        print ('No extension')

    # Save the output
    rfc.save(fc)
Harun
sumber
Apakah ada batasan ukuran file untuk Raster Split Tool? Raster saya adalah 20,8 GB! Saya terus "gagal mengeksekusi". Semua file memiliki proyeksi yang sama. Running ArcMap 10.1
Derelict
Anda selalu dapat menjalankannya dalam mode geoprocessing latar belakang 64 bit: resources.arcgis.com/en/help/main/10.1/index.html#//…
Aaron