python中grid[pos_Python有效创建密度图

[email protected],下面是一些使用罢工数据作为驱动程序的代码,而不是为每个罢工点迭代整个网格.我假设您以0.11度方格为中心,位于澳大利亚的中点,即使度的大小因纬度而有所不同！

一些可以快速参考Wikipedia的信封计算告诉我,您的40公里距离是从北到南的±4网格平方范围,而从东到西的±5网格平方范围. (在低纬度地区下降到4平方,但是…嗯！)

如前所述,这里的技巧是以直接,公式化的方式将打击位置(纬度/经度)转换为方格.找出网格的一个角的位置,从笔触中减去该位置,然后除以网格的大小-0.11度,截断,然后获得行/列索引.现在访问所有周围的正方形,直到距离变得太大为止,检查距离最多为1(2 * 2 * 4 * 5)= 81个正方形.在范围内增加平方.

结果是我最多要进行81次访问乘以1000次罢工(或者您有多少次),而不是访问100,000个网格平方乘以1000次罢工.这是显着的性能提升.

请注意,您没有描述传入的数据格式,因此我只是随机生成的数字.您将要修复该问题. 😉

#!python3

"""

Per WikiPedia (https://en.wikipedia.org/wiki/Centre_points_of_Australia)

Median point

============

The median point was calculated as the midpoint between the extremes of

latitude and longitude of the continent.

24 degrees 15 minutes south latitude, 133 degrees 25 minutes east

longitude (24°15′S 133°25′E); position on SG53-01 Henbury 1:250 000

and 5549 James 1:100 000 scale maps.

"""

MEDIAN_LAT = -(24.00 + 15.00/60.00)

MEDIAN_LON = (133 + 25.00/60.00)

"""

From the OP:

The starting grid has the overall shape of a rectangle, made up of

squares with width of 0.11 and length 0.11. The entire rectange is about

50x30. Lastly I have a shapefile which outlines the 'forecast zones' in

Australia, and if any point in the grid is outside this zone then we

omit it. So all the leftover points (insideoceanlist) are the ones in

Australia.

"""

DELTA_LAT = 0.11

DELTA_LON = 0.11

GRID_WIDTH = 50.0 # degrees

GRID_HEIGHT = 30.0 # degrees

GRID_ROWS = int(GRID_HEIGHT / DELTA_LAT) + 1

GRID_COLS = int(GRID_WIDTH / DELTA_LON) + 1

LAT_SIGN = 1.0 if MEDIAN_LAT >= 0 else -1.0

LON_SIGN = 1.0 if MEDIAN_LON >= 0 else -1.0

GRID_LOW_LAT = MEDIAN_LAT - (LAT_SIGN * GRID_HEIGHT / 2.0)

GRID_HIGH_LAT = MEDIAN_LAT + (LAT_SIGN * GRID_HEIGHT / 2.0)

GRID_MIN_LAT = min(GRID_LOW_LAT, GRID_HIGH_LAT)

GRID_MAX_LAT = max(GRID_LOW_LAT, GRID_HIGH_LAT)

GRID_LOW_LON = MEDIAN_LON - (LON_SIGN * GRID_WIDTH / 2.0)

GRID_HIGH_LON = MEDIAN_LON + (LON_SIGN * GRID_WIDTH / 2.0)

GRID_MIN_LON = min(GRID_LOW_LON, GRID_HIGH_LON)

GRID_MAX_LON = max(GRID_LOW_LON, GRID_HIGH_LON)

GRID_PROXIMITY_KM = 40.0

"""https://en.wikipedia.org/wiki/Longitude#Length_of_a_degree_of_longitude"""

_Degree_sizes_km = (

(0, 110.574, 111.320),

(15, 110.649, 107.551),

(30, 110.852, 96.486),

(45, 111.132, 78.847),

(60, 111.412, 55.800),

(75, 111.618, 28.902),

(90, 111.694, 0.000),

)

# For the Australia situation, +/- 15 degrees means that our worst

# case scenario is about 40 degrees south. At that point, a single

# degree of longitude is smallest, with a size about 80 km. That

# in turn means a 40 km distance window will span half a degree or so.

# Since grid squares a 0.11 degree across, we have to check +/- 5

# cols.

GRID_SEARCH_COLS = 5

# Latitude degrees are nice and constant-like at about 110km. That means

# a .11 degree grid square is 12km or so, making our search range +/- 4

# rows.

GRID_SEARCH_ROWS = 4

def make_grid(rows, cols):

return [[0 for col in range(cols)] for row in range(rows)]

Grid = make_grid(GRID_ROWS, GRID_COLS)

def _col_to_lon(col):

return GRID_LOW_LON + (LON_SIGN * DELTA_LON * col)

Col_to_lon = [_col_to_lon(c) for c in range(GRID_COLS)]

def _row_to_lat(row):

return GRID_LOW_LAT + (LAT_SIGN * DELTA_LAT * row)

Row_to_lat = [_row_to_lat(r) for r in range(GRID_ROWS)]

def pos_to_grid(pos):

lat, lon = pos

if lat < GRID_MIN_LAT or lat >= GRID_MAX_LAT:

print("Lat limits:", GRID_MIN_LAT, GRID_MAX_LAT)

print("Position {} is outside grid.".format(pos))

return None

if lon < GRID_MIN_LON or lon >= GRID_MAX_LON:

print("Lon limits:", GRID_MIN_LON, GRID_MAX_LON)

print("Position {} is outside grid.".format(pos))

return None

row = int((lat - GRID_LOW_LAT) / DELTA_LAT)

col = int((lon - GRID_LOW_LON) / DELTA_LON)

return (row, col)

def visit_nearby_grid_points(pos, dist_km):

row, col = pos_to_grid(pos)

# +0, +0 is not symmetric - don't increment twice

Grid[row][col] += 1

for dr in range(1, GRID_SEARCH_ROWS):

for dc in range(1, GRID_SEARCH_COLS):

misses = 0

gridpos = Row_to_lat[row+dr], Col_to_lon[col+dc]

if great_circle(pos, gridpos).meters <= dist_km:

Grid[row+dr][col+dc] += 1

else:

misses += 1

gridpos = Row_to_lat[row+dr], Col_to_lon[col-dc]

if great_circle(pos, gridpos).meters <= dist_km:

Grid[row+dr][col-dc] += 1

else:

misses += 1

gridpos = Row_to_lat[row-dr], Col_to_lon[col+dc]

if great_circle(pos, gridpos).meters <= dist_km:

Grid[row-dr][col+dc] += 1

else:

misses += 1

gridpos = Row_to_lat[row-dr], Col_to_lon[col-dc]

if great_circle(pos, gridpos).meters <= dist_km:

Grid[row-dr][col-dc] += 1

else:

misses += 1

if misses == 4:

break

def get_pos_from_line(line):

"""

FIXME: Don't know the format of your data, just random numbers

"""

import random

return (random.uniform(GRID_LOW_LAT, GRID_HIGH_LAT),

random.uniform(GRID_LOW_LON, GRID_HIGH_LON))

with open("strikes.data", "r") as strikes:

for line in strikes:

pos = get_pos_from_line(line)

visit_nearby_grid_points(pos, GRID_PROXIMITY_KM)