In [2]:
from osgeo import ogr, osr, gdal
In [3]:
# opening the geotiff file
ds = gdal.Open('G:\BTP\Satellite\Data\Test2\LE07_L1GT_147040_20050506_20170116_01_T2\LE07_L1GT_147040_20050506_20170116_01_T2_B1.TIF')
In [4]:
col, row, band = ds.RasterXSize, ds.RasterYSize, ds.RasterCount
print(col, row, band)
In [8]:
xoff, a, b, yoff, d, e = ds.GetGeoTransform()
print(xoff, a, b, yoff, d, e)
# details about the params: GDAL affine transform parameters
# xoff,yoff = left corner
# a,e = weight,height of pixels
# b,d = rotation of the image (zero if image is north up)
In [9]:
def pixel2coord(x, y):
"""Returns global coordinates from coordinates x,y of the pixel"""
xp = a * x + b * y + xoff
yp = d * x + e * y + yoff
return(xp, yp)
In [10]:
x,y = pixel2coord(col/2,row/2)
print (x, y)
In [11]:
# get the existing coordinate system
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)
In [12]:
# converting into geographic coordinate system
lonx, latx, z = transform.TransformPoint(x,y)
In [13]:
print (latx, lonx, z)
In [10]:
# rb = ds.GetRasterBand(1)
px,py = col/2,row/2 # the pixel location
pix = ds.ReadAsArray(px,py,1,1)
print pix[0][0] # pixel value
We want the name of the DISTRICT.
In [11]:
coordinates = (latx,lonx)
results = rg.search(coordinates)
print results
print type(results)
print type(results[0])
results[0]
Out[11]:
In [12]:
k = 4 # If we want k*k pixels in total from the image
for i in range(0,col,col/k):
for j in range(0,row,row/k):
# fetching the lat and lon coordinates
x,y = pixel2coord(i,j)
lonx, latx, z = transform.TransformPoint(x,y)
# fetching the name of district
coordinates = (latx,lonx)
results = rg.search(coordinates)
# The pixel value for that location
px,py = i,j
pix = ds.ReadAsArray(px,py,1,1)
pix = pix[0][0]
# printing
s = "The pixel value for the location Lat: {0:5.1f}, Long: {1:5.1f} ({2:15}) is {3:7}".format(latx,lonx,results[0]["name"],pix)
print (s)
In [13]:
g = geocoder.google([latx,lonx], method='reverse')
print type(g)
print g
print g.city
print g.state
print g.state_long
print g.country
print g.country_long
print g.address
In [14]:
results = Geocoder.reverse_geocode(latx, lonx)
print results.city
print results.country
print results.street_address
print results.administrative_area_level_1
print results.administrative_area_level_2 ## THIS GIVES THE DISTRICT !! <----------------
print results.administrative_area_level_3
In [15]:
## Converting the unicode string to ascii string
v = results.country
print type(v)
v = v.encode("ascii")
print type(v)
print v
In [16]:
k = 4 # If we want k*k pixels in total from the image
for i in range(0,col,col/k):
for j in range(0,row,row/k):
# fetching the lat and lon coordinates
x,y = pixel2coord(i,j)
lonx, latx, z = transform.TransformPoint(x,y)
# fetching the name of district
results = Geocoder.reverse_geocode(latx, lonx)
# The pixel value for that location
px,py = i,j
pix = ds.ReadAsArray(px,py,1,1)
pix = pix[0][0]
# printing
if results.country.encode('ascii') == 'India':
s = "Lat: {0:5.1f}, Long: {1:5.1f}, District: {2:12}, Pixel Val: {3:7}".format(latx,lonx,results.administrative_area_level_2,pix)
print (s)
In [17]:
import requests # To make the REST API Call
import json
In [18]:
(latx,lonx)
Out[18]:
In [19]:
url = "http://dev.virtualearth.net/REST/v1/Locations/"
point = str(latx)+","+str(lonx)
key = "Aktjg1X8bLQ_KhLQbVueYMhXDEMo7OaTweIkBvFojInYE4tVxoTp1bGKWbtU_OPJ"
response = requests.get(url+point+"?key="+key)
print(response.status_code)
In [22]:
data = response.json()
print(type(data))
In [23]:
data
Out[23]:
In [24]:
s = data["resourceSets"][0]["resources"][0]["address"]["adminDistrict2"]
s = s.encode("ascii")
s
Out[24]:
In [25]:
url = "http://dev.virtualearth.net/REST/v1/Locations/"
key = "Aktjg1X8bLQ_KhLQbVueYMhXDEMo7OaTweIkBvFojInYE4tVxoTp1bGKWbtU_OPJ"
In [33]:
k = 10 # If we want k*k pixels in total from the image
for i in range(0,col,col/k):
for j in range(0,row,row/k):
############### fetching the lat and lon coordinates #######################################
x,y = pixel2coord(i,j)
lonx, latx, z = transform.TransformPoint(x,y)
############### fetching the name of district ##############################################
point = str(latx)+","+str(lonx)
response = requests.get(url+point+"?key="+key)
data = response.json()
s = data["resourceSets"][0]["resources"][0]["address"]
if s["countryRegion"].encode("ascii") != "India":
print ("Outside Indian Territory")
continue
district = s["adminDistrict2"].encode("ascii")
############### The pixel value for that location ##########################################
px,py = i,j
pix = ds.ReadAsArray(px,py,1,1)
pix = pix[0][0]
# printing
s = "Lat: {0:5.1f}, Long: {1:5.1f}, District: {2:12}, Pixel Val: {3:7}".format(latx,lonx,district,pix)
print (s)
In [11]:
import fiona
from shapely.geometry import Point, shape
# Change this for Win7
base = "/Users/macbook/Documents/BTP/Satellite/Data/Maps/Districts/Census_2011"
fc = fiona.open(base+"/2011_Dist.shp")
In [12]:
def reverse_geocode(pt):
for feature in fc:
if shape(feature['geometry']).contains(pt):
return feature['properties']['DISTRICT']
return "NRI"
In [14]:
k = 10 # If we want k*k pixels in total from the image
for i in range(0,col,col/k):
for j in range(0,row,row/k):
############### fetching the lat and lon coordinates #######################################
x,y = pixel2coord(i,j)
lonx, latx, z = transform.TransformPoint(x,y)
############### fetching the name of district ##############################################
point = Point(lonx,latx)
district = reverse_geocode(point)
if district=="NRI":
print ("Outside Indian Territory")
continue
############### The pixel value for that location ##########################################
px,py = i,j
pix = ds.ReadAsArray(px,py,1,1)
pix = pix[0][0]
# printing
s = "Lat: {0:5.1f}, Long: {1:5.1f}, District: {2:12}, Pixel Val: {3:7}".format(latx,lonx,district,pix)
print (s)