tt.py
3.7 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
# coding=utf-8
#author: 4N
#createtime: 2021/7/15
#email: nheweijun@sina.com
from osgeo import gdal
from osgeo import gdal, gdalconst
from osgeo import ogr
import math
import requests
def tms(ytile, zoom):
n = 2.0 ** zoom
ytile = n - ytile - 1
return int(ytile)
import numpy
import cv2
from osgeo import gdal,osr
from osgeo.gdal import *
# Math functions taken from https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames #spellok
def deg2num(lat_deg, lon_deg, zoom):
lat_rad = math.radians(lat_deg)
n = 2.0 ** zoom
xtile = int((lon_deg + 180.0) / 360.0 * n)
ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
return (xtile, ytile)
def num2deg(xtile, ytile, zoom):
n = 2.0 ** zoom
lon_deg = xtile / n * 360.0 - 180.0
lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
lat_deg = math.degrees(lat_rad)
return (lat_deg, lon_deg)
def create_image(cls,image_type, pixel_array, quality):
if image_type.__eq__("image/jpeg") or image_type.__eq__("image/jpg"):
r, buf = cv2.imencode(".jpg", pixel_array, [cv2.IMWRITE_JPEG_QUALITY, quality])
image_out = buf.tobytes()
else:
height, width = pixel_array[:, :, 0].shape
four = numpy.zeros((height,width), dtype=int) + 255
four[pixel_array[:, :, 0] == 65536] = 0
r, buf = cv2.imencode(".png", numpy.dstack((pixel_array, four)))
image_out = buf.tobytes()
return image_out
if __name__ == '__main__':
lev = 8
extent = [107.78981494180618, 120.43553603935675, 18.870480519260582, 25.999868757737033]
url = "http://172.26.99.160:6060/DMap/Services/GDMap/TileServer/TMSService"
if extent[3] - extent[2] > extent[1] - extent[0]:
offset = ((extent[3] - extent[2]) - (extent[1] - extent[0])) / 2.0
out_extent = [extent[0]- offset, extent[1]+offset , extent[2], extent[3] ]
else:
offset = ((extent[1] - extent[0]) - (extent[3] - extent[2])) / 2.0
out_extent = [extent[0] , extent[1], extent[2] - offset, extent[3] + offset]
c1,r1 = deg2num(18.870480519260582,107.78981494180618,lev)
c2,r2 = deg2num(25.999868757737033, 120.43553603935675, lev)
print(r1 - r2)
print(c2 - c1)
lat1, lon1 = num2deg(c1, r1, lev)
lat2, lon2 = num2deg(c1+1, r1+1, lev)
print( [lon1, lat2, lon2, lat1])
lat3, lon3 = num2deg(c2, r2, lev)
lat4, lon4 = num2deg(c2+1, r2+1, lev)
print( [lon3, lat4, lon4, lat3])
dat = numpy.zeros(((r1 - r2+1)*256, (c2 - c1+1)*256, 3), dtype=int) + 256
offx = (extent[0]-lon1)/((lon4-lon1)/((c2 - c1+1)*256))
xnum = (extent[1]-extent[0])/((lon4-lon1)/((c2 - c1+1)*256))
print(offx)
print(xnum)
offy = (lat3 - extent[3])/((lat3-lat2)/((r1 - r2+1)*256))
ynum = (extent[3]-extent[2])/((lat3-lat2)/((r1 - r2+1)*256))
print(offy)
print(ynum)
for r in range(r2,r1+1):
for c in range(c1,c2+1):
url = "http://172.26.99.160:6060/DMap/Services/GDMap/TileServer/TMSService/{}/{}/{}.png".format(lev,c,r)
dataset:Dataset = gdal.Open(url, 0)
for b in range(1,4):
band:Band = dataset.GetRasterBand(b)
zaa = (r - r2)*256
zbb = (c - c1)*256
zcc = (r - r2+1)*256
zdd = (c - c1+1)*256
dat[(r - r2)*256:(r - r2+1)*256,(c - c1)*256 : (c - c1+1)*256, 2-(b-1)] = band.ReadAsArray(0,0,256,256)
dat2 = dat[int(offy):int(offy+ynum),int(offx):int(offx+xnum),:]
cv2.imwrite("t.jpg", dat, [cv2.IMWRITE_JPEG_QUALITY, 30])
cv2.imwrite("t2.jpg", dat2, [cv2.IMWRITE_JPEG_QUALITY, 30])