Siqi Dong's Blog

Promoting sharing of idea

0%

Posted on In
Symbols count in article: 89 Reading time ≈ 1 mins.

Baidu Map POI Industry Classification

http://lbsyun.baidu.com/index.php?title=lbscloud/poitags

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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
# -*- coding: utf-8 -*-
import json
import urllib
import urllib.request
import xlwt

def getPOIFromBaiduMap(lat_1,lon_1,lat_2,lon_2,las,keyword,excelSaveFolder):
    count = 0
    # Your BaiduMap API Key
    ak = ''
    keynew = str(urllib.parse.quote(keyword))
    push = excelSaveFolder+'\\'+keyword+'.xls'
    f = open(push, 'w')
    lat_count = int((lat_2 - lat_1) / las + 1)
    lon_count = int((lon_2 - lon_1) / las + 1)
    for lat_c in range(0, lat_count):
        lat_b1 = lat_1 + las * lat_c
        for lon_c in range(0, lon_count):
            lon_b1 = lon_1 + las * lon_c
            for i in range(0, 20):
                page_num = str(i)
                url = r'http://api.map.baidu.com/place/v2/search?query=' + keynew + '&bounds=' + str(lat_b1) + ',' + str(lon_b1) + ',' + str(lat_b1 + las) + ',' + str(lon_b1 + las) + '&page_size=20&page_num=' + str(page_num) + '&output=json&ak=' + ak
                count=count+1
                print(str(count))
                response = urllib.request.urlopen(url)
                data = json.load(response)
                try:
                    for item in data['results']:
                        jname = item['name']
                        print(jname)
                        jlat = item['location']['lat']
                        jlon = item['location']['lng']
                        jadd = item['address']
                        j_str = jname + ',' + str(jlat) + ',' + str(jlon) + ',' + jadd + '\n'
                        f.write(j_str)      
                        print(j_str)
                except:
                    print("Error")
    f.close()

# Example
# POI category
classArray=
['中餐厅',
'外国餐厅',
'小吃快餐店',
'蛋糕甜品店',
'咖啡厅',
'茶座',
'酒吧',
'星级酒店',
'快捷酒店',
'公寓式酒店',
'民宿',
'购物中心',
'百货商场',
'超市',
'便利店',
'家居建材',
'家电数码',
'商铺',
'市场',
'通讯营业厅',
'邮局',
'物流公司',
'售票处',
'洗衣店',
'图文快印店',
'照相馆',
'房产中介机构',
'公用事业',
'维修点',
'家政服务',
'殡葬服务',
'彩票销售点',
'宠物服务',
'报刊亭',
'公共厕所',
'步骑行专用道驿站',
'美容',
'美发',
'美甲',
'美体',
'公园',
'动物园',
'植物园',
'游乐园',
'博物馆',
'水族馆',
'海滨浴场',
'文物古迹',
'教堂',
'风景区',
'景点',
'寺庙',
'度假村',
'农家院',
'电影院',
'ktv',
'剧院',
'歌舞厅',
'网吧',
'游戏场所',
'洗浴按摩',
'休闲广场',
'体育场馆',
'极限运动场所',
'健身中心',
'高等院校',
'中学',
'小学',
'幼儿园',
'成人教育',
'亲子教育',
'特殊教育学校',
'留学中介机构',
'科研机构',
'培训机构',
'图书馆',
'科技馆',
'新闻出版',
'广播电视',
'艺术团体',
'美术馆',
'展览馆',
'文化宫',
'综合医院',
'专科医院',
'诊所',
'药店',
'体检机构',
'疗养院',
'急救中心',
'疾控中心',
'医疗器械',
'医疗保健',
'汽车销售',
'汽车维修',
'汽车美容',
'汽车配件',
'汽车租赁',
'汽车检测场',
'飞机场',
'火车站',
'地铁站',
'地铁线路',
'长途汽车站',
'公交车站',
'公交线路',
'港口',
'停车场',
'加油加气站',
'服务区',
'收费站',
'桥',
'充电站',
'路侧停车位',
'普通停车位',
'接送点',
'银行',
'ATM',
'信用社',
'投资理财',
'典当行',
'写字楼',
'住宅区',
'宿舍',
'内部楼栋',
'公司',
'园区',
'农林园艺',
'厂矿',
'中央机构',
'各级政府',
'行政单位',
'公检法机构',
'涉外机构',
'党派团体',
'福利机构',
'政治教育机构',
'社会团体',
'民主党派',
'居民委员会',
'高速公路出口',
'高速公路入口',
'机场出口',
'机场入口',
'车站出口',
'车站入口',
'门',
'停车场出入口',
'自行车高速出口',
'自行车高速入口',
'自行车高速出入口',
'岛屿',
'山峰',
'水系',
'省',
'省级城市',
'地级市',
'区县',
'商圈',
'乡镇',
'村庄']

# Large rectangular coordinate
lat_1=21.502705
lon_1=110.321145
lat_2=25.519951
lon_2=116.220296

# Step
las=0.2

# Get 中餐厅 POI
getPOIFromBaiduMap(lat_1, lon_1, lat_2, lon_2, las, classArray[0])

# Get all categories of POI
for i in range(len(classArray)):
    print("-------------------Collecting"+classArray[i]+"!---------------------------")
    getPOIFromBaiduMap(lat_1, lon_1, lat_2, lon_2, las, classArray[i])
    print("-------------------" + classArray[i] + "Collected-------------------------")

Posted on Edited on In
Symbols count in article: 32 Reading time ≈ 1 mins.

Data Source

http://www.gscloud.cn/

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
# -*- coding: utf-8 -*-
import time
from selenium import webdriver
from selenium.webdriver.common.keys import Keys

# chromedriver.exe path
driver=webdriver.Chrome(executable_path='XXX\\chromedriver.exe')
driver.get('http://www.gscloud.cn/accounts/login_user')
email=driver.find_element_by_xpath('//*[@id="email"]')
# Your gscloud email address
email.send_keys('XXX')
password=driver.find_element_by_xpath('//*[@id="password"]')
# Your gscloud password
password.send_keys('XXX')

captcha=driver.find_element_by_xpath('//*[@id="id_captcha_1"]')
captcha_sj=input('Please enter the code：').strip()
captcha.send_keys(captcha_sj)

dr_buttoon=driver.find_element_by_xpath('//*[@id="login-form"]/input[2]').click()
time.sleep(3)
sjzy=driver.find_element_by_link_text('数据资源').click()
time.sleep(3)
sjzy=driver.find_element_by_link_text('公开数据').click()
time.sleep(3)
sjzy=driver.find_element_by_link_text('DEM 数字高程数据').click()
time.sleep(3)
sjzy=driver.find_element_by_link_text('GDEMV2 30M 分辨率数字高程数据').click()
time.sleep(3)

# 2261 pages in all
page_num=2261
page=1
while page<=page_num:
    print('Downloading {} page'.format(page))
    # Must be 3 ~12.
    for tr_num in range(3,13):
        d_everypage='//*[@id="datasource"]/div/table/tr['+str(tr_num)+']/td[9]/div/div/p[2]'
        download=driver.find_element_by_xpath(d_everypage).click()
        # Give 20 seconds per download time
        time.sleep(20)
    page += 1
    page_sr=driver.find_element_by_xpath('//*[@id="pager1"]/div[2]/table/tr/td[7]/input')
    page_sr.clear()
    page_sr.send_keys(page)
    page_sr.send_keys(Keys.RETURN)
    time.sleep(3)

Posted on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
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
# -*- coding: utf-8 -*-
import arcpy

# The path of the DEM
demPath=""
pathList = []
elevationCell = {}
cellXSzie = arcpy.GetRasterProperties_management(demPath, "CELLSIZEX")
cellYSzie = arcpy.GetRasterProperties_management(demPath, "CELLSIZEY")
leftPosition = arcpy.GetRasterProperties_management(demPath, "LEFT")
topPosition = arcpy.GetRasterProperties_management(demPath, "TOP")
bottomPosition = arcpy.GetRasterProperties_management(demPath, "BOTTOM")
# Converting raster to array
rstArray = arcpy.RasterToNumPyArray(arcpy.Raster(demPath))
# Number of rows and columns of raster
elevationCell = {}
rows, cols = rstArray.shape
for rowNum in xrange(rows):
    for colNum in xrange(cols):
        # Set the elevation raster pixel, initialized as not flooded
        elevationCell[(rowNum, colNum)] = [rowNum, colNum, rstArray.item(rowNum, colNum), -1]

def getDirectionOfWaterFlow(rowNumber,colNumber,elevationCell):
    cell8list = []
    for rowi in range(rowNumber - 1, rowNumber + 2):
        for coli in range(colNumber - 1, colNumber + 2):
            # No out of bounds
            if (rowi < rows and coli < cols and coli >= 0 and rowi >= 0 ):
                # Exclude oneself
                if(rowi != rowNumber or coli !=colNumber):
                    # Exclude the already flooded
                    if  elevationCell[rowi, coli][3] == -1:
                        # Get elevation
                        elevationcell = elevationCell[rowi, coli]
                        cell8list.append(elevationcell)
    if len(cell8list)!=0:
        cell8list.sort(key=lambda x: (x[2]), reverse=False)
        nercelllist = cell8list
        # Get the lowest elevation grid
        mincell=nercelllist[0]
        mincell_ele = nercelllist[0][2]
        selfcell_ele = elevationCell[rowNum, colNum][2]
        # whether it can find a lower grid or not
        if mincell_ele <= selfcell_ele:
            depth = 0
            # Set to Flooded
            elevationCell[rowNum, colNumber][3] = 1
        else:
            depth = mincell_ele - selfcell_ele
            # Set to Flooded
            elevationCell[rowNumber, colNumber][3] = 1
        return mincell, depth
    else:
        return None

def getWaterFlowPath(rowNumber,colNumber,elevationCell):
    if getDirectionOfWaterFlow(rowNumber,colNumber,elevationCell)!=None:
        mincell, depth=getDirectionOfWaterFlow(rowNumber,colNumber,elevationCell)
        r=mincell[0]
        c=mincell[1]
        pathList.append(mincell)
        getWaterFlowPath(r,c,elevationCell)
    else:
        return pathList

def getWaterFlowPathRaster(pointPath,pathList,waterFlowPathRasterSaveName):
    # Set the flooded order
    for i in range(len(pathList)):
        r = pathList[i][0]
        c = pathList[i][1]
        elevationCell[(r, c)][3] = i
    for rowNum in xrange(rows):
                for colNum in xrange(cols):
                    rstArray[(rowNum, colNum)] = elevationCell[(rowNum, colNum)][3]
    result = arcpy.NumPyArrayToRaster(rstArray, arcpy.Point(leftPosition[0], bottomPosition[0]), float(cellXSzie[0]), float(cellYSzie[0]))
    result.save(waterFlowPathRasterSaveName)
    # Save the raster
    dsc = arcpy.Describe(pointPath)
    coord_sys = dsc.spatialReference
    arcpy.DefineProjection_management(waterFlowPathRasterSaveName, coord_sys)

Posted on Edited on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
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
# -*- coding: utf-8 -*-
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True

def getLocationOfPointInGrid(rasterPath,pointPath,worksapcePath):
    env.workspace = worksapcePath
    cellXSzie = arcpy.GetRasterProperties_management(rasterPath, "CELLSIZEX")
    cellYSzie = arcpy.GetRasterProperties_management(rasterPath, "CELLSIZEY")
    leftPosition = arcpy.GetRasterProperties_management(rasterPath, "LEFT")
    bottomPosition = arcpy.GetRasterProperties_management(rasterPath, "BOTTOM")
    # Converting raster to array
    rstArray = arcpy.RasterToNumPyArray(arcpy.Raster(rasterPath))
    # Number of rows and columns of raster
    rows, cols = rstArray.shape
    # Create row and column number grid
    rowcell = {}
    colcell = {}
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            rowcell[(rowNum, colNum)] = [rowNum]
            colcell[(rowNum, colNum)] = [colNum]
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            rstArray[(rowNum, colNum)] = rowcell[(rowNum, colNum)][0]
    result = arcpy.NumPyArrayToRaster(rstArray, arcpy.Point(leftPosition[0], bottomPosition[0]), float(cellXSzie[0]),float(cellYSzie[0]))
    result.save("row.tif")
    dsc = arcpy.Describe(pointPath)
    coord_sys = dsc.spatialReference
    arcpy.DefineProjection_management("row.tif", coord_sys)
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            rstArray[(rowNum, colNum)] = colcell[(rowNum, colNum)][0]
    result = arcpy.NumPyArrayToRaster(rstArray, arcpy.Point(leftPosition[0], bottomPosition[0]), float(cellXSzie[0]),float(cellYSzie[0]))
    result.save("col.tif")
    dsc = arcpy.Describe(pointPath)
    coord_sys = dsc.spatialReference
    arcpy.DefineProjection_management("col.tif", coord_sys)
    # Get the row and column of the point
    inPointFeatures = pointPath
    inRasterList = [["row.tif", "row"], ["col.tif", "col"]]
    # Execute ExtractValuesTo Points
    ExtractMultiValuesToPoints(inPointFeatures, inRasterList, "")
    cursor = arcpy.SearchCursor(inPointFeatures)
    rowN = ''
    colN = ''
    for rowcol in cursor:
        rowN = int(rowcol.row)
        colN = int(rowcol.col)
    return [rowN, colN]

Posted on Edited on In
Symbols count in article: 677 Reading time ≈ 1 mins.

The trend analysis is based on one-dimensional linear regression trend analysis, which can simulate the trend of each raster, and then integrated to reflect the spatial and temporal changes of the whole region, calculated by the following formula.

where Slope is the slope of the pixel r regression equation, the variable i is the time series number from 1 to n, n is the time span, and Mr,i is the maximum r value in the nth period.
If Slope > 0, it means that the change in r over time shows an upward trend, and the larger the Slope value, the more obvious the upward trend, and conversely, if Slope < 0, it means that the change in r over time shows a downward trend.

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
# -*- coding: utf-8 -*-
import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")

def rasterSlope(rasterFolder,resultSlopeSavePath):
	folderin = rasterFolder
	arcpy.env.workspace = folderin
	rlist = arcpy.ListRasters()
	N=len(rlist)
	i = 0
	sum1 = 0
	sum2 = 0
	sum3 = 0
	sum4 = 0
	for r in rlist:
		i += 1
		print(i)
		sum1 += i * Raster(r)
		sum2 += Raster(r)
		sum3 += i * i
		sum4 += i
		print(r)
	result = (N * sum1 - ((N + 1) * N / 2) * sum2) / (N * sum3 - sum4 * sum4)
	result.save(resultSlopeSavePath)

Posted on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# -*- coding: utf-8 -*-
import arcpy
import os

def extractSingleBand(multiBandImgPath,singleBandSavePath):
    f_list = os.listdir(multiBandImgPath)
    for file in f_list:
        if os.path.splitext(file)[1] == ".tif":
            tif = str(multiBandImgPath) + "\\" + str(file)
            arcpy.env.workspace = tif
            for raster in arcpy.ListRasters():
                #The name of band to extract
                if str(raster) == "bandName":
                    bandrasteroutputname = singleBandSavePath+"bandName" + file
                    arcpy.CopyRaster_management(raster, bandrasteroutputname, "", -999, -999)
    print  "done"

Posted on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
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
# -*- coding: utf-8 -*-
import arcpy
import os
from arcpy import env
arcpy.CheckOutExtension("Spatial")
env.overwriteOutput=1

def getExtentOfRasterOrShp(filePath):
    folder_path, file_name = os.path.split(filePath)
    arcpy.env.workspace =folder_path
    ext=file_name[-4:]
    if ext==".tif":
        extent=arcpy.Raster(filePath).extent
        return str(extent.XMin)+" "+str(extent.YMin)+" "+str(extent.XMax)+" "+str(extent.YMax)
    elif ext==".shp":
        xmin=[]
        xmax=[]
        ymin=[]
        ymax=[]
        with arcpy.da.SearchCursor(filePath, 'SHAPE@') as cursor:
            for row in cursor:
                result = row[0].extent
                xmin.append(result.XMin)
                xmax.append(result.XMax)
                ymin.append( result.YMin)
                ymax.append(result.YMax)
        xmin.sort()
        ymin.sort()
        xmax.sort(reverse=True)
        ymax.sort(reverse=True)
        return str(xmin[0])+" "+str(ymin[0])+" "+str(xmax[0])+" "+str(ymax[0])

Posted on Edited on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
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
# -*- coding: utf-8 -*-
import arcpy
import os
import sys
from arcpy.sa import *
from arcpy import env
reload(sys)
sys.setdefaultencoding('utf-8')
arcpy.CheckOutExtension("Spatial")
env.overwriteOutput=1

def getExtentOfRasterOrShp(filePath):
    folder_path, file_name = os.path.split(filePath)
    arcpy.env.workspace =folder_path
    ext=file_name[-4:]
    if ext==".tif":
        extent=arcpy.Raster(filePath).extent
        return str(extent.XMin)+" "+str(extent.YMin)+" "+str(extent.XMax)+" "+str(extent.YMax)
    elif ext==".shp":
        xmin=[]
        xmax=[]
        ymin=[]
        ymax=[]
        with arcpy.da.SearchCursor(filePath, 'SHAPE@') as cursor:
            for row in cursor:
                result = row[0].extent
                xmin.append(result.XMin)
                xmax.append(result.XMax)
                ymin.append(result.YMin)
                ymax.append(result.YMax)
        xmin.sort()
        ymin.sort()
        xmax.sort(reverse=True)
        ymax.sort(reverse=True)
        return str(xmin[0])+" "+str(ymin[0])+" "+str(xmax[0])+" "+str(ymax[0])

def MultiFieldIDWInterpolation(extentShp,shpFilePath,fieldArray,cellSize,power,searchRadius,resultSaveFolder):
    arcpy.env.workspace =resultSaveFolder
    # Get extent
    extent=getExtentOfRasterOrShp(extentShp)
    arcpy.env.extent=extent
    fields = arcpy.ListFields(shpFilePath)
    for f in fields:
        if f.name in fieldArray:
            fieldname=fieldArray[fieldArray.index(f.name)]
            # Execute IDW
            outIDW = Idw(shpFilePath, fieldname,cellSize, power, searchRadius)
            # Save the output
            outIDW.save(resultSaveFolder+"\\"+fieldname+"IDW.tif")
    print "done!"

Posted on Edited on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
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
# -*- coding: utf-8 -*-
import arcpy
import math
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")

def F_TestForTrendAnalysis(dataFolder,resultSaveFolder):
    workSpace1 = dataFolder
    arcpy.env.workspace = unicode(workSpace1, "utf8")
    ys = arcpy.ListRasters("1*", "tif")
    # Generate time series layers
    i = 0
    for y in ys:
        i = i + 1
        out = dataFolder+"\\"
        if i < 10:
            out1 = unicode(out, "utf8") + "20" + str(i) + ".tif"
        else:
            out1 = unicode(out, "utf8") + "2" + str(i) + ".tif"
        outx = Con(1, i, y)
        outx.save(out1)
    # Read the time series layer
    arcpy.env.workspace = unicode(workSpace1, "utf8")
    xs = arcpy.ListRasters("2*", "tif")
    # Calculate
    x_ = CellStatistics(xs, "MEAN", "DATA")
    y_ = CellStatistics(ys, "MEAN", "DATA")
    n = len(xs)
    xy = 0
    x2 = 0
    St = 0
    for i in range(0, len(xs)):
        St = (Raster(ys[i]) - y_) ** 2 + St
        xy = Raster(xs[i]) * Raster(ys[i]) + xy
        x2 = Raster(xs[i]) ** 2 + x2
    k = (xy - n * x_ * y_) / (x2 - n * (x_ ** 2))
    outpath1 = resultSaveFolder+"\\"
    sty1 = "Trend" + ".tif"
    out1 = outpath1 + sty1
    k.save(out1)
    # F-test
    Sr = (x2 - n * (x_ ** 2)) * (k ** 2)
    Se = St - Sr
    F = (n - 2) * Sr / Se
    outpath2 = resultSaveFolder+"\\"
    sty2 = "F" + ".tif"
    out2 = outpath2 + sty2
    F.save(out2)

Posted on Edited on In
Symbols count in article: 0 Reading time ≈ 1 mins. 
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
# -*- coding: utf-8 -*-
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("spatial")
env.overwriteOutput= 1

# get Sum
def getSumOfRasters(rasterPath,resultSavePath):
    arcpy.env.workspace = rasterPath
    rasterlist = arcpy.ListRasters("*", "TIF")
    baseraster=rasterlist[0]
    baseraster2= Con(IsNull(baseraster),0,0)
    for raster in rasterlist:
        baseraster2 = baseraster2 + raster
    baseraster2.save(resultSavePath)
    print "finish!"

# get Mean
def getMeanOfRasters(rasterPath,resultSavePath):
    arcpy.env.workspace = rasterPath
    rasterlist = arcpy.ListRasters("*", "TIF")
    baseraster=rasterlist[0]
    baseraster2= Con(IsNull(baseraster),0,0)
    for raster in rasterlist:
        baseraster2 = baseraster2 + raster
    meanraster=baseraster2/len(rasterlist)
    meanraster.save(resultSavePath)
    print "finish!"