test_code_exmaple_ttt.zip_GDAL_GPS_characterkcj_draw3u3_python

import math import utm def dd2dms(longitude, latitude): # math.modf() splits whole number and decimal into tuple # eg 53.3478 becomes (0.3478, 53) split_degx = math.modf(longitude) # the whole number [index 1] is the degrees degrees_x = int(split_degx[1]) # multiply the decimal part by 60: 0.3478 * 60 = 20.868 # split the whole number part of the total as the minutes: 20 # abs() absoulte value - no negative minutes_x = abs(int(math.modf(split_degx[0] * 60)[1])) # multiply the decimal part of the split above by 60 to get the seconds # 0.868 x 60 = 52.08, round excess decimal places to 2 places # abs() absoulte value - no negative seconds_x = abs(round(math.modf(split_degx[0] * 60)[0] * 60, 2)) # repeat for latitude split_degy = math.modf(latitude) degrees_y = int(split_degy[1]) minutes_y = abs(int(math.modf(split_degy[0] * 60)[1])) seconds_y = abs(round(math.modf(split_degy[0] * 60)[0] * 60, 2)) # account for E/W & N/S if degrees_x < 0: EorW = "W" else: EorW = "E" if degrees_y < 0: NorS = "S" else: NorS = "N" # abs() remove negative from degrees, was only needed for if-else above print("\t" + str(abs(degrees_x)) + u"\u00b0 " + str(minutes_x) + "' " + str(seconds_x) + "\" " + EorW) print("\t" + str(abs(degrees_y)) + u"\u00b0 " + str(minutes_y) + "' " + str(seconds_y) + "\" " + NorS) # some coords of cities coords = [["Dublin", -6.2597, 53.3478], ["Paris", 2.3508, 48.8567], ["Sydney", 151.2094, -33.8650]] # test dd2dms() for city, x, y in coords: print(city + ":", dd2dms(x, y)) # x is longitude, y is latitude tif_filename = "/data/MSC_080428180432_09355_20161270_1G_PS/MSC_080428180432_09355_20161270_PS.tif" # reference : https://gis.stackexchange.com/questions/60371/gdal-python-how-do-i-get-coordinate-system-name-from-spatialreference from osgeo import gdal,osr ds=gdal.Open(tif_filename) prj=ds.GetProjection() print(prj) srs=osr.SpatialReference(wkt=prj) if srs.IsProjected: print(srs.GetAttrValue('projcs')) print(srs.GetAttrValue('geogcs')) from osgeo import osr, gdal # get the existing coordinate system ds = gdal.Open(tif_filename) 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"]]""" wgs84_wkt = prj 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) transform = osr.CoordinateTransformation(old_cs,old_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) # print(latlong) # reference : http://monkut.webfactional.com/blog/archive/2012/5/2/understanding-raster-basic-gis-concepts-and-the-python-gdal-library/ # RasterIO : https://geohackweek.github.io/raster/04-workingwithrasters/ from osgeo import gdal,ogr,osr def GetExtent(gt,cols,rows): ''' Return list of corner coordinates from a geotransform @type gt: C{tuple/list} @param gt: geotransform @type cols: C{int} @param cols: number of columns in the dataset @type rows: C{int} @param rows: number of rows in the dataset @rtype: C{[float,...,float]} @return: coordinates of each corner ''' ext=[] xarr=[0,cols] yarr=[0,rows] for px in xarr: for py in yarr: x=gt[0]+(px*gt[1])+(py*gt[2]) y=gt[3]+(px*gt[4])+(py*gt[5]) ext.append([x,y]) print(x,y) yarr.reverse() return ext def ReprojectCoords(coords,src_srs,tgt_srs): ''' Reproject a list of x,y coordinates. @type geom: C{tuple/list} @param geom: List of [[x,y],...[x,y]] coordinates @type src_srs: C{osr.SpatialReference} @param src_srs: OSR SpatialReference object @type tgt_srs: C{osr.SpatialReference} @param tgt_srs: OSR SpatialReference object @rtype: C{tuple/list} @return: List of transformed [[x,y],...[x,y]] coordinates ''' trans_coords=[] transform = osr.CoordinateTransformation( src_srs, tgt_srs) for x,y in coords: x,y,z = transform.TransformPoint(x,y) trans_coords.append([x,y]) return trans_coords raster=tif_filename ds=gdal.Open(raster) gt=ds.GetGeoTransform() cols = ds.RasterXSize rows = ds.RasterYSize ext=GetExtent(gt,cols,rows) src_srs=osr.SpatialReference() src_srs.ImportFromWkt(ds.GetProjection()) #tgt_srs=osr.SpatialReference() #tgt_srs.ImportFromEPSG(4326) tgt_srs = src_srs.CloneGeogCS() geo_ext=ReprojectCoords(ext,src_srs,tgt_srs) print(geo_ext) # Get UTM import subprocess imagefile = tif_filename p= subprocess.Popen(["gdalinfo", "%s"%imagefile], stdout=subprocess.PIPE) out,err= p.communicate() outs = out.decode(encoding="utf-8") print(type(outs)) print(outs) out = outs ul= out[out.find("Upper Left")+15:out.find("Upper Left")+38] print(ul) lr= out[out.find("Lower Right")+15:out.find("Lower Right")+38] print(lr) center = out[out.find("Center")+15:out.find("Center")+38] print(center) utm_zone = out[out.find("UTM zone ")+9:out.find("UTM zone ")+12] print(utm_zone) easting, northing = center.split(', ') # print(type(utm_x), type(utm_y)) import re utm_zone_str = re.findall(r'\d+|[^\s]', utm_zone) print(utm_zone_str) # UTM to DD import utm # utm.to_latlon(EASTING, NORTHING, ZONE NUMBER, ZONE LETTER). print('utm to latlon', utm.to_latlon(float(easting), float(northing), int(utm_zone_str[0]), utm_zone_str[1])) latlon = utm.to_latlon(float(easting), float(northing), int(utm_zone_str[0]), utm_zone_str[1]) latlon = list(latlon) print(latlon) # DD to DMS print(dd2dms(latlon[1], latlon[0])) # x is longitude, y is latitude