tạo trường cao độ/chiều cao gdal gump numpy python

Question

Tôi muốn tạo một số độ cao/độ cao rasters bằng cách sử dụng python, gdal và gumpy. Tôi đang mắc kẹt trên NumPy (. Và có lẽ python và GDAL)

Trong NumPy, tôi đã cố gắng như sau:

>>> a= numpy.linspace(4,1,4, endpoint=True) 
>>> b= numpy.vstack(a) 
>>> c=numpy.repeat(b,4,axis=1) 
>>> c 
array([[ 4., 4., 4., 4.], 
     [ 3., 3., 3., 3.], 
     [ 2., 2., 2., 2.], 
     [ 1., 1., 1., 1.]]) #This is the elevation data I want 

từ OSGeo nhập khẩu GDAL từ gdalconst nhập khẩu *

>>> format = "Terragen" 
>>> driver = gdal.GetDriverByName(format)  
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE = 1',                'MAXUSERPIXELVALUE= 4']) 
>>> raster = numpy.zeros((4,4), dtype=numpy.float32) #this is where I'm messing up 
>>> dst_ds.GetRasterBand(1).WriteArray(raster) # gives the null elevation data I asked for in (raster) 
0 
>>> dst_ds = None 

Tôi nghĩ rằng tôi đang thiếu một cái gì đó đơn giản và mong được lời khuyên của bạn.

Cảm ơn,

Chris

(tiếp theo sau)

• terragendataset.cpp, v 1,2 *
  • dự án: Terragen (tm) TER driver
  • Mục đích: Tài liệu đọc cho Terragen TER
  • Tác giả: Ray Gardener, Daylon Graphics TNHH *
  • phần của mô-đun này có nguồn gốc từ trình điều khiển GDAL bởi
  • Frank Warmerdam, xem http://www.gdal.org

lời xin lỗi của tôi trước để Ray Gardener và Frank Warmerdam.

Terragen định dạng ghi chú:

Đối viết: scal = gridpost khoảng cách tính bằng mét hv_px = hv_m/scal span_px = span_m/scal bù đắp = thấy TerragenDataset :: write_header() quy mô = thấy TerragenDataset: : write_header() hv vật lý = (hv_px - bù đắp) * 65.536,0/quy mô

Chúng tôi nói với người gọi rằng:

Elevations are Int16 when reading, 
    and Float32 when writing. We need logical 
    elevations when writing so that we can 
    encode them with as much precision as possible 
    when going down to physical 16-bit ints. 
    Implementing band::SetScale/SetOffset won't work because 
    it requires callers to know format write details. 
    So we've added two Create() options that let the 
    caller tell us the span's logical extent, and with 
    those two values we can convert to physical pixels. 

      ds::SetGeoTransform() lets us establish the 
     size of ground pixels. 
    ds::SetProjection() lets us establish what 
     units ground measures are in (also needed 
     to calc the size of ground pixels). 
    band::SetUnitType() tells us what units 
     the given Float32 elevations are in. 

này nói với tôi rằng trước khi WriteArray trên của tôi (somearray) tôi phải đặt cả Biến đổi địa và setProjection và SetUnitType để mọi thứ hoạt (khả năng)

Từ Tutorial GDAL API: Python nhập khẩu OSR nhập numpy

dst_ds.SetGeoTransform([ 444720, 30, 0, 3751320, 0, -30 ]) 

srs = osr.SpatialReference() 
srs.SetUTM(11, 1) 
srs.SetWellKnownGeogCS('NAD27') 
dst_ds.SetProjection(srs.ExportToWkt()) 

raster = numpy.zeros((512, 512), dtype=numpy.uint8) #we've seen this before 
dst_ds.GetRasterBand(1).WriteArray(raster) 

# Once we're done, close properly the dataset 
dst_ds = None 

Tôi xin lỗi vì đã tạo bài đăng quá dài và thú nhận. Nếu và khi tôi nhận được quyền này, nó sẽ là tốt đẹp để có tất cả các ghi chú ở một nơi (bài dài).

Các Confession:

Tôi đã từng chụp ảnh (jpeg), chuyển đổi nó vào một GeoTIFF, và nhập khẩu nó như gạch vào một cơ sở dữ liệu PostGIS. Tôi hiện đang cố gắng tạo raster độ cao để treo lên hình ảnh. Điều này có vẻ ngớ ngẩn, nhưng có một nghệ sĩ mà tôi muốn tôn vinh, trong khi tại cùng một thời gian không xúc phạm những người đã làm việc một cách tỉ mỉ để tạo và cải thiện những công cụ tuyệt vời này.

Nghệ sĩ là người Bỉ nên mét sẽ theo thứ tự. Cô làm việc ở hạ Manhattan, NY, UTM 18. Bây giờ một số SetGeoTransform hợp lý. Hình ảnh được đề cập ở trên là w = 3649/h = 2736, tôi sẽ phải suy nghĩ về điều này.

Một thử:

>>> format = "Terragen" 
>>> driver = gdal.GetDriverByName(format) 
>>> dst_ds = driver.Create('test_3.ter', 4,4,1, gdal.GDT_Float32, ['MINUSERPIXELVALUE=1','MAXUSERPIXELVALUE-4']) 
>>> type(dst_ds) 
<class 'osgeo.gdal.Dataset'> 
>>> import osr 
>>> dst_ds.SetGeoTransform([582495, 1, 0.5, 4512717, 0.5, -1]) #x-res 0.5, y_res 0.5 a guess 
0 
>>> type(dst_ds) 
<class 'osgeo.gdal.Dataset'> 
>>> srs = osr.SpatialReference() 
>>> srs.SetUTM(18,1) 
0 
>>> srs.SetWellKnownGeogCS('WGS84') 
0 
>>> dst_ds.SetProjection(srs.ExportToWkt()) 
0 
>>> raster = c.astype(numpy.float32) 
>>> dst_ds.GetRasterBand(1).WriteArray(raster) 
0 
>>> dst_ds = None 
>>> from osgeo import gdalinfo 
>>> gdalinfo.main(['foo', 'test_3.ter']) 
Driver: Terragen/Terragen heightfield 
Files: test_3.ter 
Size is 4, 4 
Coordinate System is: 
LOCAL_CS["Terragen world space", 
    UNIT["m",1]] 
Origin = (0.000000000000000,0.000000000000000) 
Pixel Size = (1.000000000000000,1.000000000000000) 
Metadata: 
    AREA_OR_POINT=Point 
Corner Coordinates: 
Upper Left ( 0.0000000, 0.0000000) 
Lower Left ( 0.0000000, 4.0000000) 
Upper Right ( 4.0000000, 0.0000000) 
Lower Right ( 4.0000000, 4.0000000) 
Center  ( 2.0000000, 2.0000000) 
Band 1 Block=4x1 Type=Int16, ColorInterp=Undefined 
    Unit Type: m 
Offset: 2, Scale:7.62939453125e-05 
0 

Rõ ràng tiến gần hơn nhưng rõ ràng nếu SetUTM (18,1) được vớt lên. Đây có phải là 4x4 trong số Sông Hudson hoặc Local_CS (hệ thống tọa độ) không? Một sự im lặng là gì?

More Reading

// Terragen files aren't really georeferenced, but 
// we should get the projection's linear units so 
// that we can scale elevations correctly. 

// Increase the heightscale until the physical span 
// fits within a 16-bit range. The smaller the logical span, 
// the more necessary this becomes. 

4x4 mét là một khoảng logic khá nhỏ.

Vì vậy, có lẽ điều này là tốt như nó được. Các SetGeoTransform được các đơn vị phải, thiết lập quy mô và bạn có không gian thế giới Terragen của bạn.

Tư tưởng cuối cùng: Tôi không lập trình, nhưng đến mức tôi có thể theo dõi. Điều đó nói rằng, tôi kinda băn khoăn đầu tiên và sau đó nếu những gì dữ liệu trông giống như trong Terragen nhỏ của tôi Thế giới Space (nhờ sau để http://www.gis.usu.edu/~chrisg/python/2009/ tuần 4):

>>> fn = 'test_3.ter' 
>>> ds = gdal.Open(fn, GA_ReadOnly) 
>>> band = ds.GetRasterBand(1) 
>>> data = band.ReadAsArray(0,0,1,1) 
>>> data 
array([[26214]], dtype=int16) 
>>> data = band.ReadAsArray(0,0,4,4) 
>>> data 
array([[ 26214, 26214, 26214, 26214], 
     [ 13107, 13107, 13107, 13107], 
     [  0,  0,  0,  0], 
     [-13107, -13107, -13107, -13107]], dtype=int16) 
>>> 

Vì vậy, đây là hài lòng. Tôi tưởng tượng sự khác biệt giữa trên được sử dụng numpy c và kết quả này đi đến các hành động áp dụng Float16 qua này rất nhỏ khoảng logic.

Và vào 'hf2':

>>> src_ds = gdal.Open('test_3.ter') 
>>> dst_ds = driver.CreateCopy('test_6.hf2', src_ds, 0) 
>>> dst_ds.SetGeoTransform([582495,1,0.5,4512717,0.5,-1]) 
0 
>>> srs = osr.SpatialReference() 
>>> srs.SetUTM(18,1) 
0 
>>> srs.SetWellKnownGeogCS('WGS84') 
0 
>>> dst_ds.SetProjection(srs.ExportToWkt()) 
0 
>>> dst_ds = None 
>>> src_ds = None 
>>> gdalinfo.main(['foo','test_6.hf2']) 
Driver: HF2/HF2/HFZ heightfield raster 
Files: test_6.hf2 
    test_6.hf2.aux.xml 
Size is 4, 4 
Coordinate System is: 
PROJCS["UTM Zone 18, Northern Hemisphere", 
    GEOGCS["WGS 84", 
     DATUM["WGS_1984", 
      SPHEROID["WGS 84",6378137,298.257223563, 
      AUTHORITY["EPSG","7030"]], 
     TOWGS84[0,0,0,0,0,0,0], 
     AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich",0, 
     AUTHORITY["EPSG","8901"]], 
    UNIT["degree",0.0174532925199433, 
     AUTHORITY["EPSG","9108"]], 
    AUTHORITY["EPSG","4326"]], 
PROJECTION["Transverse_Mercator"], 
PARAMETER["latitude_of_origin",0], 
PARAMETER["central_meridian",-75], 
PARAMETER["scale_factor",0.9996], 
PARAMETER["false_easting",500000], 
PARAMETER["false_northing",0], 
UNIT["Meter",1]] 
Origin = (0.000000000000000,0.000000000000000) 
Pixel Size = (1.000000000000000,1.000000000000000) 
Metadata: 
VERTICAL_PRECISION=1.000000 
Corner Coordinates: 
Upper Left ( 0.0000000, 0.0000000) (79d29'19.48"W, 0d 0' 0.01"N) 
Lower Left ( 0.0000000, 4.0000000) (79d29'19.48"W, 0d 0' 0.13"N) 
Upper Right ( 4.0000000, 0.0000000) (79d29'19.35"W, 0d 0' 0.01"N) 
Lower Right ( 4.0000000, 4.0000000) (79d29'19.35"W, 0d 0' 0.13"N) 
Center  ( 2.0000000, 2.0000000) (79d29'19.41"W, 0d 0' 0.06"N) 
Band 1 Block=256x1 Type=Float32, ColorInterp=Undefined 
Unit Type: m 
0 
>>> 

Gần hoàn toàn hài lòng, dù có vẻ như tôi đang ở La Concordia Peru. Vì vậy, tôi có để tìm hiểu cách nói - giống như miền Bắc hơn, như New York North. Có SetUTM lấy một phần tử thứ ba cho thấy 'cách xa' bắc hay nam. Dường như tôi đã chạy qua một biểu đồ UTM ngày hôm qua đã có khu nhãn thư, một cái gì đó giống như C tại đường xích đạo và có thể T hoặc S cho khu vực New York. Tôi thực sự nghĩ rằng SetGeoTransform về cơ bản là thiết lập phía trên bên trái phía bắc và phía đông và điều này đã ảnh hưởng đến phần 'xa về phía bắc/nam' của SetUTM như thế nào. Tắt đến gdal-dev.

Và sau đó vẫn còn:

Paddington Bear đến Peru vì anh ấy có vé. Tôi đến đó vì tôi nói đó là nơi tôi muốn. Terragen, làm việc như nó, đã cho tôi không gian pixel của tôi. Các cuộc gọi tiếp theo để srs hành động theo hf2 (SetUTM), nhưng phía đông và phía bắc được thành lập theo Terragen để UTM 18 đã được thiết lập, nhưng trong một hộp giới hạn tại đường xích đạo. Đủ tốt.

gdal_translate đưa tôi đến New York. Tôi đang ở trong cửa sổ nên một dòng lệnh.và kết quả;

C:\Program Files\GDAL>gdalinfo c:/python27/test_9.hf2 
Driver: HF2/HF2/HFZ heightfield raster 
Files: c:/python27/test_9.hf2 
Size is 4, 4 
Coordinate System is: 
PROJCS["UTM Zone 18, Northern Hemisphere", 
GEOGCS["NAD83", 
    DATUM["North_American_Datum_1983", 
     SPHEROID["GRS 1980",6378137,298.257222101, 
      AUTHORITY["EPSG","7019"]], 
     TOWGS84[0,0,0,0,0,0,0], 
     AUTHORITY["EPSG","6269"]], 
    PRIMEM["Greenwich",0, 
     AUTHORITY["EPSG","8901"]], 
    UNIT["degree",0.0174532925199433, 
     AUTHORITY["EPSG","9122"]], 
    AUTHORITY["EPSG","4269"]], 
PROJECTION["Transverse_Mercator"], 
PARAMETER["latitude_of_origin",0], 
PARAMETER["central_meridian",-75], 
PARAMETER["scale_factor",0.9996], 
PARAMETER["false_easting",500000], 
PARAMETER["false_northing",0], 
UNIT["Meter",1]] 
Origin = (583862.000000000000000,4510893.000000000000000) 
Pixel Size = (-1.000000000000000,-1.000000000000000) 
Metadata: 
VERTICAL_PRECISION=0.010000 
Corner Coordinates: 
Upper Left ( 583862.000, 4510893.000) (74d 0'24.04"W, 40d44'40.97"N) 
Lower Left ( 583862.000, 4510889.000) (74d 0'24.04"W, 40d44'40.84"N) 
Upper Right ( 583858.000, 4510893.000) (74d 0'24.21"W, 40d44'40.97"N) 
Lower Right ( 583858.000, 4510889.000) (74d 0'24.21"W, 40d44'40.84"N) 
Center  ( 583860.000, 4510891.000) (74d 0'24.13"W, 40d44'40.91"N) 
Band 1 Block=256x1 Type=Float32, ColorInterp=Undefined 
Unit Type: m 

C:\Program Files\GDAL> 

Vì vậy, chúng tôi trở lại NY. Có lẽ cách tốt hơn để tiếp cận tất cả điều này. Tôi phải có mục tiêu đã chấp nhận Tạo vì tôi đã đăng tải/ứng biến một tập dữ liệu từ numpy. Tôi cần phải xem xét các định dạng khác cho phép tạo. Độ cao trong GeoTiff là một khả năng (tôi nghĩ vậy.)

Cảm ơn tất cả các nhận xét, đề xuất và nudges nhẹ nhàng đối với việc đọc thích hợp. Tạo bản đồ trong python thật thú vị!

Chris

Answer 1

Bạn không quá xa. Vấn đề duy nhất của bạn là vấn đề cú pháp với các tùy chọn tạo GDAL. Thay thế:

['MINUSERPIXELVALUE = 1','MAXUSERPIXELVALUE= 4'] 

với không có khoảng trống trước hoặc sau các cặp khóa/giá trị:

['MINUSERPIXELVALUE=1', 'MAXUSERPIXELVALUE=4'] 

Kiểm tra type(dst_ds) và nó phải <class 'osgeo.gdal.Dataset'> hơn <type 'NoneType'> cho một im lặng thất bại như là trường hợp cho ở trên.

---

Cập nhật Theo mặc định, warnings or exceptions are not shown. Bạn có thể bật chúng qua gdal.UseExceptions() để xem nội dung nào đang đánh dấu bên dưới, ví dụ:

>>> from osgeo import gdal 
>>> gdal.UseExceptions() 
>>> driver = gdal.GetDriverByName('Terragen') 
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE = 1', 'MAXUSERPIXELVALUE= 4']) 
Warning 6: Driver Terragen does not support MINUSERPIXELVALUE creation option 
>>> dst_ds = driver.Create('test.ter', 4,4,1,gdal.GDT_Float32, ['MINUSERPIXELVALUE=1', 'MAXUSERPIXELVALUE=4']) 
>>> type(dst_ds) 
<class 'osgeo.gdal.Dataset'>