I have my own python application that I'm thinking of making. For that purpose, it is better to have a map displayed. It's like a study for that. I hope it helps someone.
There seems to be a way to use openstreetmap and basemap module, folium, shp file and GIS software, but I didn't seem to find it difficult or suitable for me. After all, when I looked up the google static map API, it was easier than I expected and the aerial photography suits my purpose, so I decided to use google map. It seems that the Geographical Survey Institute map can be handled in the same way, so it's like ...
Regarding the conversion between latitude / longitude and the web Mercator projection, TRAIL NOTE The author's explanation was easy to understand. There are some things that I haven't fully understood about coordinate transformation, etc., but that point may become apparent as the development progresses further.
--I checked the operation with Win10 / 32bit Anaconda from Spyder and a little earlier version of raspian Jessie. --Operates and displays on both PC screens and low resolution screens for raspberry pi --The tiles displayed when online can be saved as a file and can be displayed offline. --For the time being, the display is openCV. This just adopted the method that I understand at the moment. It connects, cuts, displays, saves, and reopens tiles. At present, there is nothing special about image processing. Regarding opencv installation --For windows here --For raspberrypi, I tried to use opencv3 + python3.x by referring to here, but I gave up because it was insufficient for 8G disk. It seems that the usage of urllib is different, so some code is rewritten and it works with opncv2 + python2 --For the time being, only keyboard operation ---: Shrink, +: Enlarge, a, s, d, f: Move to west, north, south, east respectively --Space, backspace: Switch map type (google hybrid → google aerial photograph → google roadmap → Geographical Survey Institute standard map →…) ――It takes a little preparation to actually operate it. Dig a directory for saving tiles
import os
os.mkdir("maptiles")
for d in range(22):
os.mkdir("maptiles/z%02d"%d)
Looking at it like this, it may be Akan to put a long slapstick on it ... I will be careful in future posts.
AMAP.py
# -*- coding: utf-8 -*-
import numpy as np
import cv2
import urllib
import os
#WIN_W=480
#WIN_H=260
WIN_W=1280
WIN_H=960
# gg*:googleMap, other:cyberjapandata:Geospatial Information Authority of Japan
map_type_name = ["gghybrid","ggsatellite","ggroadmap", "std","ort_old10","gazo1","seamlessphoto"]
TILE_W = [640,640,640,256,256,256,256]
TILE_H = [640,640,640,256,256,256,256]
min_zoom=[ 0, 0, 0, 2,10,10, 2]
max_zoom=[21,21,21,18,17,17,18]
fmt=["png","png","png", "png", "png","jpg","jpg"]
#You can see various maps by changing the table definition
##ort:2007-, airphoto:2004-, gazo4:1988-1990, gazo3:1984-1986, gazo2:1979-1983, gazo1:1974-1978
##ort_old10:1961-1964, ort_USA10:1945-1950
#map_type_name = ["ort","airphoto","gazo4","gazo3","gazo2","gazo1","ort_old10","ort_USA10"]
#TILE_W = [256,256,256,256,256,256,256,256]
#TILE_H = [256,256,256,256,256,256,256,256]
#min_zoom=[14, 5,10,10,10,10,10,10]
#max_zoom=[18,18,17,17,17,17,17,17]
#fmt=["jpg","png","jpg","jpg","jpg","jpg","png","png"]
#Tokyo Station
HOME_LON=139.767052
HOME_LAT= 35.681167
HOME_ZOOM=18
TILES_DIR="maptiles/"
max_pixels= [256*2**zm for zm in range(22)]
#A dictionary that stores references to open tiles. Keyed by map type, zoom, index x, index y
opened_tiles={}
white_tiles={}
#lon:longitude, lat:latitude
def ll2pix(lon, lat, zoom):
pix_x=2**(zoom+7)*(lon/180+1)
pix_y=2**(zoom+7)*(-np.arctanh(np.sin(np.pi/180*lat))/np.pi+1)
return pix_x,pix_y
def pix2ll(x,y,zoom):
lon=180*(x/(2**(zoom+7))-1)
lat=180/np.pi*(np.arcsin(np.tanh(-np.pi/(2**(zoom+7))*y+np.pi)))
return lon, lat
#From longitude / latitude dx,dy Returns the longitude / latitude when moving pixels
def new_ll(lon_cur,lat_cur,zm, dx,dy):
x,y=ll2pix(lon_cur,lat_cur,zm)
return pix2ll(x+dx,y+dy,zm)
def dddmm2f(dddmm_mmmm):
#12345.6789 ->123 degrees 45.6789 minutes->123 degrees.(45.6789/60)
ddd=int(dddmm_mmmm)//100
mm_mmmm=dddmm_mmmm-ddd*100
return ddd+mm_mmmm/60
#Concatenate tiles to create an image larger than the display window, and finally cut it to the window size and return it
def load_win_img(mtype, lon,lat,zm):
cx,cy=ll2pix(lon,lat,zm)
win_left=int(cx-WIN_W/2)
win_top=int(cy-WIN_H/2)
x_nth=win_left//TILE_W[mtype]
y_nth=win_top//TILE_H[mtype]
left_offset = win_left%TILE_W[mtype]
top_offset = win_top%TILE_H[mtype]
vcon_list=[]
tot_height=0
tot_height += TILE_H[mtype]-top_offset
j=0
while True:
hcon_list=[]
tot_width=0
tot_width += TILE_W[mtype]-left_offset
i=0
while True:
img_tmp=open_tile_img(mtype, x_nth+i,y_nth+j,zm)
hcon_list.append(img_tmp) #
if tot_width >= WIN_W:
break
tot_width += TILE_W[mtype]
i+=1
hcon_img=cv2.hconcat(hcon_list)
vcon_list.append(hcon_img)
if tot_height >= WIN_H:
break
tot_height += TILE_H[mtype]
j+=1
convined_img=cv2.vconcat(vcon_list)
return convined_img[top_offset:top_offset+WIN_H, left_offset:left_offset+WIN_W, :]
def tile_file_name(mtype, x_nth,y_nth,zm):
# x_nth=x//TILE_W[mtype]
# y_nth=y//TILE_H[mtype]
return TILES_DIR+"z%02d/%s_z%02d_%dx%d_%07d_%07d"%(zm,map_type_name[mtype],zm,TILE_W[mtype],TILE_H[mtype],x_nth,y_nth)+"."+fmt[mtype]
#Open tile
#A point that has never been opened->Try to download / save. If it fails, it returns a whitewashed image. If successful, register the point / zoom / map type in the dictionary and mark it as opened thereafter.
#Saved to a file but not open->Open normally. Also register in the dictionary
#Already open->Returns the image registered in the dictionary
def open_tile_img(mtype, x_nth,y_nth,zm):
if (mtype, zm,x_nth,y_nth) in opened_tiles:
print("opened_tiles(%d,%d,%d,%d)"%(mtype, zm,x_nth,y_nth))
return opened_tiles[(mtype, zm,x_nth,y_nth)]
fname=tile_file_name(mtype, x_nth,y_nth,zm)
if os.path.exists(fname):
print("opening tile(%d,%d,%d,%d)"%(mtype,zm,x_nth,y_nth) +" -> "+fname)
else:
c_lon,c_lat=pix2ll((x_nth+0.5)*TILE_W[mtype],(y_nth+0.5)*TILE_H[mtype],zm)
if map_type_name[mtype][0:2]=="gg":
url="http://maps.google.com/maps/api/staticmap?"
url+="¢er=%.08f,%08f&zoom=%d&size=%dx%d&maptype=%s" % \
(c_lat,c_lon,zm,TILE_W[mtype],TILE_H[mtype],map_type_name[mtype][2:])
#maptype
#roadmap Normal map. Default value for maptype parameter
#satellite aerial photography
#terrain Physical terrain map image showing terrain and vegetation
#hybrid aerial photograph + normal map. Layered main roads and place names on top of aerial photographs
else:
url="http://cyberjapandata.gsi.go.jp/xyz/%s/%d/%d/%d.%s"%(map_type_name[mtype],zm,x_nth,y_nth,fmt[mtype])
print("Downloading... ")
print(url)
print(" -> "+fname)
try:
urllib.request.urlretrieve(url,fname) #python3
# urllib.urlretrieve(url,fname) #python2
except Exception as e:
#If the tile cannot be obtained, the image filled with white is returned.
print(e)
print("Download faild -> blank")
if (TILE_W[mtype],TILE_H[mtype]) in white_tiles:
return white_tiles[(TILE_W[mtype],TILE_H[mtype])]
else:
white=np.zeros([TILE_H[mtype],TILE_W[mtype],3],dtype=np.uint8)
white[:,:,:]=255
white_tiles[(TILE_W[mtype],TILE_H[mtype])]=white
return white
opened_tiles[(mtype, zm,x_nth,y_nth)]=cv2.imread(fname)
return opened_tiles[(mtype, zm,x_nth,y_nth)]
if __name__ == '__main__':
map_type=0
c_lon=HOME_LON
c_lat=HOME_LAT
zoom=HOME_ZOOM
cv2.namedWindow("Ackerman's Map", cv2.WINDOW_AUTOSIZE)
map_type_bak = -1
c_lon_bak = -1
c_lat_bak = -1
zoom_bak = -1
#mainloop
while (True):
if map_type_bak != map_type or c_lon_bak != c_lon or c_lat_bak != c_lat or zoom_bak != zoom:
win_img=load_win_img(map_type, c_lon,c_lat,zoom)
cv2.imshow("Ackerman's Map", win_img)
map_type_bak = map_type
c_lon_bak = c_lon
c_lat_bak = c_lat
zoom_bak = zoom
k=cv2.waitKey(0) & 0xff #Polling without setting the waiting time to 0 when adding the GPS function
print("pressed:"+str(k))
if k == ord('+'):
if zoom<max_zoom[map_type]:
zoom += 1
elif k == ord('-'):
if zoom>min_zoom[map_type]:
zoom -= 1
elif k == ord('a'):
c_lon,c_lat=new_ll(c_lon,c_lat,zoom,-WIN_W/4,0)
elif k == ord('s'):
c_lon,c_lat=new_ll(c_lon,c_lat,zoom,0,-WIN_H/4)
elif k == ord('d'):
c_lon,c_lat=new_ll(c_lon,c_lat,zoom,0,+WIN_H/4)
elif k == ord('f'):
c_lon,c_lat=new_ll(c_lon,c_lat,zoom,+WIN_W/4,0)
elif k == 32:
#space key
map_type = (map_type+1)%len(map_type_name)
if zoom > max_zoom[map_type]:
zoom=max_zoom[map_type]
if zoom < min_zoom[map_type]:
zoom=min_zoom[map_type]
elif k == 8:
#Backspace
map_type = (map_type-1)%len(map_type_name)
if zoom > max_zoom[map_type]:
zoom=max_zoom[map_type]
if zoom < min_zoom[map_type]:
zoom=min_zoom[map_type]
elif k == ord("q"):
break
cv2.destroyAllWindows()
For a while, I traced the example of the introductory book of raspberry pi and modified the program from time to time, and tried a short code while reading the introductory book of python3, but I wonder why I wrote it from scratch for the first time. I used to program mainly in C language more than 15 years ago, but compared to that situation, it's a list, a dictionary, an inclusive notation, and heaven. What is the usefulness of associative arrays (≒ dictionaries) in other languages? I was thinking, but I feel that I was able to use it effectively here.
Also, the Geographical Survey Institute map displayed as a side note is quite interesting. "Did you live here 50 years ago?" "Wow, this mysterious building has existed for more than 30 years ?!" In the code
#Tokyo Station
HOME_LON=139.767052
HOME_LAT= 35.681167
It might be interesting to take a walk while switching the map by changing the area around to the latitude and longitude of your residence.
I've already written about connecting and displaying my current location, and I've made something like a degraded version of FoxtrotGPS, but it was published in a magazine. I haven't posted it for the time being because the code is quite packed. It's like deciphering the NMEA format, but it was easier than I expected. I was also considering a method like communicating with gpsd, but it seems easier than that.
Recommended Posts