ऊंचाई/ऊंचाई क्षेत्र GDAL numpy अजगर

Question

मैं कुछ ऊंचाई/heightfield अजगर, GDAL और numpy का उपयोग कर rasters बनाना चाहेंगे बनाने। मैं numpy पर (। अजगर और GDAL और शायद)

numpy में अटक कर रहा हूँ, मैं प्रयास कर रहा है निम्नलिखित:

>>> 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 

osgeo आयात GDAL से gdalconst आयात *

>>> 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 

से

मैं समझ मैं कुछ सरल याद कर रहा हूँ और आपकी सलाह के लिए तत्पर हैं।

धन्यवाद,

क्रिस

(बाद में जारी रखा)

• terragendataset.cpp, वी 1,2 *
  • परियोजना: Terragen (टीएम) TER ड्राइवर
  • उद्देश्य: Terragen TER के लिए रीडर दस्तावेजों
  • लेखक: रे माली, Daylon ग्राफिक्स लिमिटेड *
  • इस मॉड्यूल
  • फ्रैंक Warmerdam द्वारा GDAL चालकों से प्राप्त के भाग, http://www.gdal.org

मेरे रे माली और फ्रैंक के लिए अग्रिम में क्षमा याचना देखना Warmerdam।

Terragen प्रारूप नोट:

लेखन के लिए: scal = मीटर में gridpost दूरी hv_px = hv_m/scal span_px = span_m/scal ऑफसेट = TerragenDataset देख :: write_header() पैमाने = TerragenDataset देखें: : write_header() शारीरिक एचवी = (hv_px - ऑफसेट) * 65,536.0/पैमाने

हम कॉल करने बताते हैं कि:

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. 

यह मुझसे कहता है कि मेरे ऊपर WriteArray (somearray) मैं सेट दोनों के लिए जियोट्रांस्फ़ॉर्म और SetProjection और चीजों को काम करने के लिए (संभावित) SetUnitType

GDAL एपीआई ट्यूटोरियल से है करने से पहले: अजगर आयात OSR आयात वाले अत्यंत लंबे पद और एक बयान बनाने के लिए

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 

मेरे क्षमायाचना Numpy। और जब मैं यह अधिकार मिलता है, यह एक ही स्थान पर (लंबी पोस्ट) में सभी नोट करने के लिए अच्छा होगा।

बयान:

मैं पहले एक तस्वीर (JPEG) ले लिया है, एक GeoTiff करने के लिए इसे परिवर्तित, और एक PostGIS डेटाबेस में टाइल्स के रूप में यह आयात किया। अब मैं तस्वीर को ढेर करने के लिए ऊंचाई रास्टर बनाने का प्रयास कर रहा हूं। यह शायद मूर्खतापूर्ण प्रतीत होता है, लेकिन एक कलाकार है जिसे मैं सम्मान करना चाहता हूं, जबकि पर उन लोगों को अपमानित नहीं किया जिन्होंने इन महान उपकरणों को बनाने और सुधारने के लिए दृढ़ता से काम किया है।

कलाकार बेल्जियम है तो मीटर क्रम में होगा। वह निचले मैनहट्टन, एनवाई, यूटीएम 18 में काम करती है। अब कुछ उचित SetGeoTransform। उपर्युक्त तस्वीर w = 3649/h = 2736 है, मुझे इसे कुछ सोचना होगा।

एक और कोशिश:

>>> 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 

जाहिर करीब लेकिन यह स्पष्ट नहीं है, तो SetUTM (18,1) उठाया गया था हो रही। क्या यह हडसन नदी या स्थानीय_CS (समन्वय प्रणाली) में 4x4 है? एक मूक असफल क्या है?

अधिक पढ़ना

// 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 मीटर एक बहुत छोटा सा तार्किक अवधि है।

तो, शायद यह उतना ही अच्छा है जितना इसे प्राप्त होता है। SetGeoTransform इकाइयों को सही बनाता है, पैमाने निर्धारित करता है और आपके पास टेरेगन वर्ल्ड स्पेस है।

अंतिम विचार: मैं प्रोग्राम नहीं करता हूं, लेकिन कुछ हद तक मैं साथ चल सकता हूं। जैसा कि कहा गया, मैं थोड़े पहले अगर और फिर क्या डेटा मेरी छोटी Terragen वर्ल्ड स्पेस (http://www.gis.usu.edu/~chrisg/python/2009/ सप्ताह 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) 
>>> 

तो यह संतुष्टिदायक है। मैं उपरोक्त प्रयुक्त numpy c के बीच अंतर की कल्पना करता हूं और यह परिणाम फ्लोट 16 को इस बहुत छोटे लॉजिकल अवधि में लागू करने के कार्यों पर जाता है।

और '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 
>>> 

लगभग पूरी तरह से संतोषजनक है, हालांकि ऐसा लगता है कि मैं ला Concordia पेरू में हूँ। तो मेरे पास है कि यह कहने के लिए कि कैसे अधिक उत्तर की तरह, न्यूयॉर्क उत्तर की तरह। क्या SetUTM उत्तर या दक्षिण 'कितना दूर' का सुझाव देने वाला तीसरा तत्व लेता है। ऐसा लगता है कि कल मैंने एक यूटीएम चार्ट में भाग लिया था जिसमें पत्र लेबल क्षेत्र थे, भूमध्य रेखा पर सी जैसे कुछ और शायद न्यूयॉर्क क्षेत्र के लिए टी या एस।

मैंने वास्तव में सोचा था कि SetGeoTransform अनिवार्य रूप से शीर्ष बाएं उत्तर और पूर्व की स्थापना कर रहा था और यह SetUTM के 'उत्तर/दक्षिण' भाग को प्रभावित कर रहा था। Gdal-dev के लिए बंद।

और बाद में अभी भी:

पेडिंगटन बीयर पेरू के लिए चला गया, क्योंकि वह एक टिकट था। मैं वहां गया क्योंकि मैंने कहा था कि वह वही था जहां मैं बनना चाहता था। टेरेगन, जैसा काम करता है, ने मुझे अपनी पिक्सेल स्पेस दी। एसआरएस के बाद की कॉल ने एचएफ 2 (सेटयूटीएम) पर काम किया, लेकिन पूर्व और उत्तर की स्थापना टेरेगेन के तहत की गई थी ताकि यूटीएम 18 सेट हो, लेकिन भूमध्य रेखा पर एक बाउंडिंग बॉक्स में। काफी है।

gdal_translate मुझे न्यूयॉर्क ले गया। मैं खिड़कियों में एक कमांड लाइन में हूँ।और नतीजा;

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> 

तो, हम NY में वापस आ गए हैं। इस सब से संपर्क करने के लिए शायद बेहतर तरीके हैं। I को एक ऐसा लक्ष्य होना था जो स्वीकार किया गया था क्योंकि मैं डेटासेट को पोस्टसेट/सुधार कर रहा था। मुझे अन्य प्रारूपों को देखने की ज़रूरत है जो बनाने की अनुमति देते हैं। जियोटीफ में ऊंचाई एक संभावना है (मुझे लगता है।)

उचित पढ़ने की दिशा में सभी टिप्पणियों, सुझावों और सौम्य nudges के लिए धन्यवाद। पाइथन में मानचित्र बनाना मजेदार है!

क्रिस

Answer 1

आप बहुत दूर नहीं हैं। आपकी एकमात्र समस्या जीडीएएल निर्माण विकल्पों के साथ वाक्यविन्यास समस्या है। बदलें:

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

साथ

कोई रिक्त स्थान से पहले या कुंजी/मान जोड़े के बाद:

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

चेक type(dst_ds) और यह <class 'osgeo.gdal.Dataset'> बजाय <type 'NoneType'> होना चाहिए एक मूक के रूप में भी यही स्थिति थी असफल के लिए ऊपर।

---

अद्यतन डिफ़ॉल्ट रूप से, warnings or exceptions are not shown। आप gdal.UseExceptions() के माध्यम से नीचे क्या टिक रहा है यह देखने के लिए सक्षम कर सकते हैं, उदाहरण:

>>> 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'>