2010-05-27 21 views
30

पायथन में जीडीएएल का उपयोग करके, आप GeoTIFF फ़ाइल का अक्षांश और देशांतर कैसे प्राप्त करते हैं?एक GeoTIFF फ़ाइल से अक्षांश और देशांतर प्राप्त करें

GeoTIFF किसी भी समन्वय जानकारी को स्टोर करने के लिए प्रकट नहीं होता है। इसके बजाय, वे एक्सवाई मूल निर्देशांक स्टोर करते हैं। हालांकि, एक्सवाई निर्देशांक शीर्ष बाएं कोने और नीचे बाएं कोने के अक्षांश और देशांतर प्रदान नहीं करते हैं।

ऐसा प्रतीत होता है कि मुझे इस समस्या को हल करने के लिए कुछ गणित करने की आवश्यकता होगी, लेकिन मेरे पास कहां से शुरू करना है इस पर कोई सुराग नहीं है।

ऐसा करने के लिए किस प्रक्रिया की आवश्यकता है?

मुझे पता है कि GetGeoTransform() विधि इसके लिए महत्वपूर्ण है, हालांकि, मुझे नहीं पता कि वहां से इसके साथ क्या करना है।

उत्तर

59

अपने GeoTiff के कोनों के निर्देशांक प्राप्त करने निम्न करें:

from osgeo import gdal 
ds = gdal.Open('path/to/file') 
width = ds.RasterXSize 
height = ds.RasterYSize 
gt = ds.GetGeoTransform() 
minx = gt[0] 
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2] 
maxy = gt[3] 

हालांकि, इन अक्षांश/देशांतर प्रारूप में नहीं हो सकता है। जस्टिन ने नोट किया, आपका जियोटिफ़ किसी प्रकार की समन्वय प्रणाली के साथ संग्रहीत किया जाएगा। आप क्या प्रणाली यह है समन्वय नहीं जानते हैं, तो आप gdalinfo चलाकर पता कर सकते हैं:

gdalinfo ~/somedir/somefile.tif 

कौन सा आउटपुट:

Driver: GTiff/GeoTIFF 
Size is 512, 512 
Coordinate System is: 
PROJCS["NAD27/UTM zone 11N", 
    GEOGCS["NAD27", 
     DATUM["North_American_Datum_1927", 
      SPHEROID["Clarke 1866",6378206.4,294.978698213901]], 
     PRIMEM["Greenwich",0], 
     UNIT["degree",0.0174532925199433]], 
    PROJECTION["Transverse_Mercator"], 
    PARAMETER["latitude_of_origin",0], 
    PARAMETER["central_meridian",-117], 
    PARAMETER["scale_factor",0.9996], 
    PARAMETER["false_easting",500000], 
    PARAMETER["false_northing",0], 
    UNIT["metre",1]] 
Origin = (440720.000000,3751320.000000) 
Pixel Size = (60.000000,-60.000000) 
Corner Coordinates: 
Upper Left ( 440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N) 
Lower Left ( 440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N) 
Upper Right ( 471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N) 
Lower Right ( 471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N) 
Center  ( 456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N) 
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray 

यह आउटपुट आप सभी की जरूरत हो सकती है। यदि आप इस प्रोग्राम को पाइथन में करना चाहते हैं, तो इस तरह, आपको वही जानकारी मिलती है।

यदि समन्वय प्रणाली PROJCS है तो ऊपर दिए गए उदाहरण की तरह आप एक अनुमानित समन्वय प्रणाली से निपट रहे हैं। एक अनुमानित समन्वयक प्रणाली representation of the spheroidal earth's surface, but flattened and distorted onto a plane है। यदि आप अक्षांश और देशांतर चाहते हैं, तो आपको निर्देशांक को भौगोलिक समन्वय प्रणाली में परिवर्तित करना होगा जो आप चाहते हैं।

अफसोस की बात है, पृथ्वी के विभिन्न गोलाकार मॉडल के आधार पर सभी अक्षांश/देशांतर जोड़े बराबर नहीं बनाए जाते हैं। इस उदाहरण में, मैं WGS84 में परिवर्तित कर रहा हूं, भौगोलिक समन्वय प्रणाली जीपीएस में अनुकूल है और सभी लोकप्रिय वेब मैपिंग साइटों द्वारा उपयोग की जाती है। समन्वय प्रणाली एक अच्छी तरह से परिभाषित स्ट्रिंग द्वारा परिभाषित किया गया है। उनमें से एक सूची spatial ref से उपलब्ध है, उदाहरण के लिए देखें WGS84

from osgeo import osr, gdal 

# get the existing coordinate system 
ds = gdal.Open('path/to/file') 
old_cs= osr.SpatialReference() 
old_cs.ImportFromWkt(ds.GetProjectionRef()) 

# create the new coordinate system 
wgs84_wkt = """ 
GEOGCS["WGS 84", 
    DATUM["WGS_1984", 
     SPHEROID["WGS 84",6378137,298.257223563, 
      AUTHORITY["EPSG","7030"]], 
     AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich",0, 
     AUTHORITY["EPSG","8901"]], 
    UNIT["degree",0.01745329251994328, 
     AUTHORITY["EPSG","9122"]], 
    AUTHORITY["EPSG","4326"]]""" 
new_cs = osr.SpatialReference() 
new_cs .ImportFromWkt(wgs84_wkt) 

# create a transform object to convert between coordinate systems 
transform = osr.CoordinateTransformation(old_cs,new_cs) 

#get the point to transform, pixel (0,0) in this case 
width = ds.RasterXSize 
height = ds.RasterYSize 
gt = ds.GetGeoTransform() 
minx = gt[0] 
miny = gt[3] + width*gt[4] + height*gt[5] 

#get the coordinates in lat long 
latlong = transform.TransformPoint(x,y) 

उम्मीद है कि यह वही करेगा जो आप चाहते हैं।

+1

वाह!बस मुझे जो कुछ चाहिए वह सब करता है। धन्यवाद! :) – Phanto

+1

अंतिम पंक्ति में, एक्स और वाई परिभाषित नहीं हैं – blaylockbk

10

अगर यह एक पूरा जवाब है मैं नहीं जानता, लेकिन this site का कहना है:

x/y नक्शा आयाम पूरब की ओर गमन और उत्तर की ओर कहा जाता है। एक भौगोलिक समन्वय प्रणाली में डेटासेट के लिए इन्हें देशांतर और अक्षांश बनाएगा। अनुमानित समन्वय प्रणाली के लिए वे आम तौर पर अनुमानित समन्वय प्रणाली में पूर्वोत्तर और उत्तरदायी होंगे। पूर्वनिर्धारित छवियों के लिए पूर्व और उत्तर केवल प्रत्येक पिक्सेल के पिक्सेल/लाइन ऑफसेट होंगे (जैसा कि एकता geotransform द्वारा निहित)।

इसलिए वे वास्तव में देशांतर और अक्षांश हो सकते हैं।

+0

मुझे विश्वास नहीं है कि मैंने उस पृष्ठ पर यह नहीं देखा था। धन्यवाद! (यदि मेरी रैंक अधिक थी, तो मैं आपकी पोस्ट को अपनाना चाहूंगा) – Phanto

संबंधित मुद्दे