أكثر

تقسيم البيانات النقطية إلى أجزاء أصغر باستخدام GDAL؟

تقسيم البيانات النقطية إلى أجزاء أصغر باستخدام GDAL؟


لديّ نقطية (USGS DEM في الواقع) وأحتاج إلى تقسيمها إلى أجزاء أصغر مثل الصورة أدناه. تم تحقيق ذلك في ArcGIS 10.0 باستخدام أداة Split Raster. أرغب في طريقة البرمجيات الحرة والمفتوحة المصدر للقيام بذلك. لقد نظرت إلى GDAL ، وأعتقد أنها ستفعل ذلك (بطريقة ما باستخدام gdal_translate) ، لكن لا يمكنني العثور على أي شيء. في النهاية ، أود أن أكون قادرًا على أخذ البيانات النقطية وأقول إلى أي مدى أريد تقسيمها إلى (4KM من 4KM قطع).


سيعمل gdal_translate باستخدام خياري -srcwin أو -projwin.

-srcwin xoff yoff xsize ysize: تحديد إطار فرعي من الصورة المصدر للنسخ بناءً على موقع البكسل / الخط.

-projwin ulx uly lrx lry: تحديد إطار فرعي من الصورة المصدر للنسخ (مثل -srcwin) ولكن مع الزوايا الواردة في الإحداثيات المحددة جغرافيًا.

ستحتاج إلى تحديد مواقع البكسل / الخط أو إحداثيات الزاوية ثم تكرار القيم باستخدام gdal_translate. سيعمل شيء مثل python السريع والقذر أدناه إذا كان استخدام قيم البكسل و -srcwin مناسبًا لك ، فسيكون هناك المزيد من العمل لفرز الإحداثيات.

استيراد نظام التشغيل ، gdal من gdalconst import * العرض = 512 الارتفاع = 512 البلاط الحجم = 64 بالنسبة لـ i في النطاق (0 ، العرض ، حجم البلاط): لـ j في النطاق (0 ، الارتفاع ، حجم البلاط): gdaltranString = "gdal_translate -of GTIFF -srcwin "+ str (i) +"، "+ str (j) +"، "+ str (البلاط) +"، " + str (البلاط) +" utm.tif utm _ "+ str (i) +" _ " + str (j) + ". tif" os.system (gdaltranString)

يقرأ الحل الذي أقدمه ، استنادًا إلى الحل منwwnick ، ​​الأبعاد النقطية من الملف نفسه ، ويغطي الصورة بأكملها بجعل مربعات الحافة أصغر حجمًا إذا لزم الأمر:

استيراد نظام التشغيل ، sys من osgeo استيراد gdal dset = gdal.Open (sys.argv [1]) العرض = dset.RasterXSize height = dset.RasterYS حجم عرض الطباعة ، 'x' ، ارتفاع البلاط الحجم = 5000 لـ i في النطاق (0 ، العرض ، حجم البلاط): بالنسبة لـ j في النطاق (0 ، الارتفاع ، حجم البلاط): w = min (i + البلاط ، العرض) - ih = min (j + البلاط ، الارتفاع) - j gdaltranString = "gdal_translate -of GTIFF -srcwin" + str (i) + "،" + str (j) + "،" + str (w) + "،"  + str (h) + "" + sys.argv [1] + "" + sys.argv [ 2] + "_" + str (i) + "_" + str (j) + ". tif" os.system (gdaltranString)

يوجد نص برمجي بيثون مجمّع خصيصًا لتجميع البيانات النقطية ، gdal_retile:

gdal_retile.py [-v] [-co NAME = VALUE] * [-of out_format] [-ps pixelWidth pixelHeight] [-overlap val_in_pixel] [-ot {Byte / Int16 / UInt16 / UInt32 / Int32 / Float32 / Float64 / CInt16 / CInt32 / CFloat32 / CFloat64}] '[-tileIndex tileIndexName [-tileIndexField tileIndexFieldName]] [-csv fileName [-csvDelim delimiter]] [-s_srs srs_def] [-pyramidOnly] [-r {near / cubic cubic lanczos}] -عدد مستويات المستويات [-useDirForEachRow] -targetDir TileDirectory input_files

على سبيل المثال:

gdal_retile.py -ps 512 512 -targetDir C: example dir some_dem.tif


يمكنك استخدام r.tile من GRASS GIS. ينشئ r.tile خريطة نقطية منفصلة لكل بلاطة بأسماء خريطة مرقمة بناءً على بادئة يحددها المستخدم. يمكن تحديد عرض البلاط (الأعمدة) وارتفاع البلاط (الصفوف).

باستخدام Python API للجلسة العشبية ، هناك حاجة فقط إلى بضعة أسطر من كود Python لاستدعاء وظيفة r.tile من الخارج ، أي لكتابة نص مستقل. باستخدام r.external و r.external.out أيضًا ، لا يحدث تكرار للبيانات أثناء خطوة معالجة GRASS GIS.

كود مزيف:

  1. تهيئة جلسة العشب
  2. تحديد تنسيق الإخراج مع r.external.out
  3. ربط ملف الإدخال مع r.external
  4. تشغيل r.tile الذي يولد المربعات بالتنسيق المحدد أعلاه
  5. قم بإلغاء ربط r.external.out
  6. جلسة العشب المغلقة

لـaron الذي سأل:

آمل أن أجد نسخة gdalwarp من إجابة @ wwnick التي تستخدم الخيار -Multi لتحسين العمليات متعددة النواة ومتعددة مؤشرات الترابط

إخلاء طفيف

هذا يستخدمغدالوارب، على الرغم من أنني لست مقتنعًا تمامًا بأنه سيكون هناك مكاسب كبيرة في الأداء. حتى الآن كنت ملتزمًا بإدخال / إخراج - تشغيل هذا البرنامج النصي على خطوط نقطية كبيرة وتقطيعه إلى العديد من الأجزاء الأصغر لا يبدو مكثفًا لوحدة المعالجة المركزية ، لذلك أفترض أن عنق الزجاجة يكمن في الكتابة على القرص. إذا كنت تخطط لإعادة إسقاط المربعات في نفس الوقت أو شيء مشابه ، فقد يتغير هذا. هناك نصائح ضبط هنا. لم تسفر مسرحية قصيرة عن أي تحسن بالنسبة لي ، ولم يكن يبدو أن وحدة المعالجة المركزية هي العامل المحدد.

إخلاء المسؤولية جانبًا ، إليك النص الذي سيتم استخدامهغدالواربلتقسيم البيانات النقطية إلى عدة بلاطات أصغر. قد يكون هناك بعض الخسارة بسبب تقسيم الأرضية ولكن يمكن التعامل مع ذلك باختيار عدد البلاط الذي تريده. سيكون ذلكن + 1أيننهو الرقم الذي تقسم عليه للحصول علىعرض_الجدولوارتفاع البلاطالمتغيرات.

استيراد عملية فرعية استيراد gdal استيراد sys def gdalwarp (* args): إرجاع subprocess.check_call (['gdalwarp'] + list (args)) src_path = sys.argv [1] ds = gdal.Open (src_path) حاول: out_base = sys .argv [2] باستثناء IndexError: out_base = '/ tmp / test_' gt = ds.GetGeoTransform () width_px = ds.RasterXSize height_px = ds.RasterYSize # احصل على أسلاك للزاوية اليسرى السفلية xmin = int (gt [0]) xmax = int (gt [0] + (gt [1] * width_px)) # احصل على أسلاك للزاوية اليمنى العليا إذا كانت gt ​​[5]> 0: ymin = int (gt [3] - (gt [5] * height_px)) else: ymin = int (gt [3] + (gt [5] * height_px)) ymax = int (gt [3]) # تقسيم الارتفاع والعرض إلى أربعة - أي أن هذا سينتج 25 بلاطة بلاط_عرض = (xmax - xmin) // 4 Tile_height = (ymax - ymin) // 4 لـ x في النطاق (xmin، xmax، Tile_width): لـ y في النطاق (ymin، ymax، tile_height): gdalwarp ('- te'، str (x)، str (y)، str (x + tile_width)، str (y + tile_height)، '-multi'، '-wo'، 'NUM_THREADS = ALL_CPUS'، '-wm'، '500'، src_path، out_base + '{} _ {}. tif'.format (x، y))

استيراد gdal import os def create_tiles (tile_size، input_filename، in_path = "/ media / Data / avinash / input /"، out_path = "/ home / nesac / PycharmProjects / GIS_Data_Automation / data / البلاط /"): ds = gdal.Open ( in_path + input_filename) band = ds.GetRasterBand (1) output_filename = 'البلاط_' xsize ، ysize = (band.XSize ، band.YSize) tile_size_x ، Tile_size_y = Tile_size tile_list = {} complete_x = xsize // tile_size_x // complete_x Tile_size_y متبقي_x = xsize٪ tile_size_x متبقية_y = ysize٪ tile_size_y # للجزء A لـ i في النطاق (complete_x): بالنسبة لـ j في النطاق (complete_y): Xmin = i * tile_size_x Xmax = i * tile_size_x + Tile_size_x - 1 Ysmin_y * Ymax = j * tile_size_y + tile_size_y - 1 # إنشاء التصحيح هنا com_string = "gdal_translate -of GTIFF -srcwin" + str (Xmin) + "،" + str (Ymin) + "،" + str (tile_size_x) + "، "+ str (tile_size_y) +" "+  str (in_path) + str (input_filename) +" "+ str (out_path) + str (output_filename) + str (Xmin) +" _ "+ str (Ymin) +". tif "os.system (com_string ) x_residue_count = 1 y_residue_count = 1 # للجزء B لـ j في النطاق (complete_y): Xmin = tile_size_x * complete_x Xmax = Tile_size_x * complete_x + Residue_x - 1 Ymin = j * tile_size_y Ymax = j * tile_size_y + 1_s الإنشاء هنا com_string = "gdal_translate -of GTIFF -srcwin" + str (Xmin) + "،" + str (Ymin) + "،" + str (المتبقية_x) + "،" + str (tile_size_y) + "" +  str (in_path) + str (input_filename) + "" + str (out_path) + str (output_filename) + str (Xmin) + "_" + str (Ymin) + ".tif" os.system (com_string) # للجزء C بالنسبة لـ i في النطاق (complete_x): Xmin = i * tile_size_x Xmax = i * tile_size_x + Tile_size_x - 1 Ymin = tile_size_y * complete_y Ymax = tile_size_y * complete_y + البقايا_ص - 1 com_string = "gdal_translate -of GTIFF-strsr ) + "،" + str (Ymin) + "،" + str (tile_size_x) + "،" + str (المتبقية_y) + "" +  str (in_path) + str (input_filename) + "" + str (out_path) + str (output_filename) + str (Xmin) + "_" + str (Ymin) + ".tif" os.system (com_strin ز) # للجزء د Xmin = complete_x * tile_size_x Ymin = complete_y * tile_size_y com_string = "gdal_translate -of GTIFF -srcwin" + str (Xmin) + "،" + str (Ymin) + "،" + str (المتبقية_x) + "،" + str (المتبقية_y) + "" +  str (in_path) + str (input_filename) + "" + str (out_path) + str (output_filename) + str (Xmin) + "_" + str (Ymin) + ".tif" os.system (com_string)

شاهد الفيديو: طريقة وصل عدة خلايا في خلية واحدة في الإكسل