<jupyter_start><jupyter_text>**Requirements/Dependecies**<jupyter_code>!pip install rasterio
import rasterio
import numpy as np
from rasterio.plot import show, adjust_band
import numpy as np
from rasterio.warp import transform
from rasterio.crs import CRS
from math import sin, cos, sqrt, atan2, radians
import math
import matplotlib.pyplot as plt<jupyter_output>Collecting rasterio
[?25l  Downloading https://files.pythonhosted.org/packages/e1/bf/d3c5e7df3828db144a6797269bf3aec31db96c20f13e75b93179eb059955/rasterio-1.2.3-cp37-cp37m-manylinux1_x86_64.whl (19.1MB)
[K     |████████████████████████████████| 19.1MB 220kB/s 
[?25hCollecting affine
  Downloading https://files.pythonhosted.org/packages/ac/a6/1a39a1ede71210e3ddaf623982b06ecfc5c5c03741ae659073159184cd3e/affine-2.3.0-py2.py3-none-any.whl
Requirement already satisfied: certifi in /usr/local/lib/python3.7/dist-packages (from rasterio) (2020.12.5)
Collecting cligj>=0.5
  Downloading https://files.pythonhosted.org/packages/42/1e/947eadf10d6804bf276eb8a038bd5307996dceaaa41cfd21b7a15ec62f5d/cligj-0.7.1-py3-none-any.whl
Requirement already satisfied: click<8,>=4.0 in /usr/local/lib/python3.7/dist-packages (from rasterio) (7.1.2)
Collecting snuggs>=1.4.1
  Downloading https://files.pythonhosted.org/packages/cc/0e/d27d6e806d6c0d1a2cfdc5d1f088e42339a0a54a09c3343f7f81ec8947ea/snuggs-1.4.7-py3-non[...]<jupyter_text>**Load Image**<jupyter_code>dataset = rasterio.open('/content/TQ28ne_DTM_2m.tif')<jupyter_output><empty_output><jupyter_text>**Coordinate refernece System**<jupyter_code>dataset.crs<jupyter_output><empty_output><jupyter_text>**Image bounds**<jupyter_code>dataset.bounds<jupyter_output><empty_output><jupyter_text>**Affine Transform Info**<jupyter_code>dataset.transform<jupyter_output><empty_output><jupyter_text>**Channels**<jupyter_code>dataset.indexes<jupyter_output><empty_output><jupyter_text>**Read Channel**<jupyter_code>img = dataset.read(1)
imgdata = np.array([adjust_band(dataset.read(1)) for i in (3,2,1)])
show(imgdata*3)<jupyter_output>Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).<jupyter_text>**Pixel to Image location in Image CRS**<jupyter_code>x, y = dataset.xy(1, 1)
print(x,y)
x, y = dataset.transform * (1000, 1000)
print(x,y)
row, column = dataset.index(527000.0, 188000.0)
print(row,column)<jupyter_output>1000 1000<jupyter_text>**Image CRS to WSG84**<jupyter_code>img_pts = ([(dataset.transform * (x, x)) for x in np.arange(1, 1000, 1)])

WGS84 = CRS.from_epsg(4326)
epsg_27700 = CRS.from_epsg(27700)

def get_coords(img_pts):
  results = []
  for pos in img_pts:
    lat_long = transform(epsg_27700, WGS84,xs=[pos[0]], ys=[pos[1]])
    results.append([lat_long[1][0],lat_long[0][0]])
  return results

lat_longs = get_coords(img_pts)
print(lat_longs)
new_crs = CRS.from_epsg(4326)
new_coo = transform(dataset.crs, new_crs,xs=[527000.0], ys=[188000.0])
new_coo
new_coo = transform(new_crs,dataset.crs,xs=[-0.19677678491595763], ys=[51.594968736823255])
new_coo<jupyter_output><empty_output><jupyter_text>**Get Elevation with Image CRS**<jupyter_code>coords = ((525200,189800),)

elevation = "/content/TQ28ne_DTM_2m.tif"
with rasterio.open(elevation) as src:
    vals = src.sample(coords)
    for val in vals:
        print(val[0])<jupyter_output>75.087494<jupyter_text>**Elevation from Lat Longs**<jupyter_code>WGS84 = CRS.from_epsg(4326)
epsg_27700 = CRS.from_epsg(27700)

def get_elevation(lat_lon,image_path):
  results = []
  #read the data
  src = rasterio.open(image_path)

  for i,ele in enumerate(lat_lon,0):
    lat = ele[0]
    lon = ele[1]
    #Conversion to 
    new_coo = transform(WGS84,epsg_27700,xs=[lon], ys=[lat])
    vals = src.sample(((new_coo[0][0],new_coo[1][0]),))
    for val in vals:
      results.append({"lat":lat,"lon":lon,"elevation":val[0]})
  return results

test_results = get_elevation(lat_longs,'/content/TQ28ne_DTM_2m.tif')
elevation_measures = [child["elevation"] for child in test_results]
elevation_measures<jupyter_output><empty_output><jupyter_text>**BASIC STAT INFO**<jupyter_code>mean_elev=round((sum(elevation_measures)/len(elevation_measures)),3)
min_elev=min(elevation_measures)
max_elev=max(elevation_measures)
R = 6373.0

lat1 = radians(lat_longs[0][0])
lon1 = radians(lat_longs[0][1])
lat2 = radians(lat_longs[-1][0])
lon2 = radians(lat_longs[-1][1])

dlon = lon2 - lon1
dlat = lat2 - lat1

a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))

distance = R * c

print("Result:", distance)
print("Should be:", 278.546, "km")<jupyter_output>Result: 2.8190371472128537
Should be: 278.546 km<jupyter_text>**Distance Measure**<jupyter_code>def haversine(lat1,lon1,lat2,lon2):
    lat1_rad=math.radians(lat1)
    lat2_rad=math.radians(lat2)
    lon1_rad=math.radians(lon1)
    lon2_rad=math.radians(lon2)
    delta_lat=lat2_rad-lat1_rad
    delta_lon=lon2_rad-lon1_rad
    a=math.sqrt((math.sin(delta_lat/2))**2+math.cos(lat1_rad)*math.cos(lat2_rad)*(math.sin(delta_lon/2))**2)
    d=2*6371000*math.asin(a)
    return d

d_list=[]
for j in range(len(lat_longs)):
    lat_p=lat_longs[j][0]
    lon_p=lat_longs[j][1]
    dp=haversine(lat_longs[0][0],lat_longs[0][1],lat_p,lon_p)/1000 #km
    d_list.append(dp)
d_list_rev=d_list[::-1]<jupyter_output><empty_output><jupyter_text>**Elevation Profile**<jupyter_code>base_reg=0
plt.figure(figsize=(10,4))
plt.plot(d_list_rev,elevation_measures)
plt.plot([0,distance],[min_elev,min_elev],'--g',label='min: '+str(min_elev)+' m')
plt.plot([0,distance],[max_elev,max_elev],'--r',label='max: '+str(max_elev)+' m')
plt.plot([0,distance],[mean_elev,mean_elev],'--y',label='ave: '+str(mean_elev)+' m')
plt.fill_between(d_list_rev,elevation_measures,base_reg,alpha=0.1)
plt.text(d_list_rev[0],elevation_measures[0],"P1")
plt.text(d_list_rev[-1],elevation_measures[-1],"P2")
plt.xlabel("Distance(km)")
plt.ylabel("Elevation(m)")
plt.grid()
plt.legend(fontsize='small')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1644451849383539, 'fraction_numerical': 0.08730423192269243, 'mean_word_length': 7.268595041322314, 'pattern_counts': {'":': 3, '<': 53, '<?xml version=': 0, '>': 55, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '80501', 'license': 'MIT', 'n_tokens_mistral': 2917, 'n_tokens_neox': 2603, 'n_words': 330, 'path': 'visualisations/altitude_profile.ipynb', 'repo_name': 'xesaad/DSGApril2021CityMaas', 'size': '6002'}