#python #scikit-learn #gis #kernel-density
#python #scikit-learn #гис #плотность ядра
Вопрос:
Используя плотность ядра sklearn.neighbors, я получаю значения плотности, которые намного меньше, чем я ожидаю. Оценки плотности имеют значения, которые составляют примерно 1/200 от того, что я ожидал.
Я рассмотрел оценку плотности ядра sklearn для распределения видов, и я пошел по этому пути преобразования моих входных данных широты / долготы в радианы и использования метрики расстояния haversine, но я получал странные результаты.
Я много думал об этом, и вот параметры, которые имеют для меня наибольший смысл.
- пропускная способность = 1,0 единицы отображения
- метрика = «Евклидово» для моего 2D-пространства
- kernel = «Gaussian» для плотности вероятности
Вопросы
- Кажется ли это разумным подходом к этой проблеме?
- Почему значения плотности намного меньше, чем я ожидал?
Вот моя функция и параметры, которые я ей передаю.
import rasterio
from rasterio.crs import CRS
from sklearn.neighbors import KernelDensity
import numpy as np
def kernel_density_lat_lon(positions, bandwidth, metric, kernel,
cell_size, extent, output_raster, multiplier=None):
# Set the bounds of the output raster based on the extent
x_min = extent[0]
x_max = extent[1]
y_min = extent[2]
y_max = extent[3]
# Create arrays, based on cell_size and bounds
# These arrays hold x locations and y locations for each pixel in the output raster
x = np.arange(x_min, x_max, cell_size)
y = np.arange(y_min, y_max, cell_size)
# Create a meshgrid, which has cells whose values are the (x,y) location at each cell
xx, yy = np.meshgrid(x, y)
# Pair the x locations with y locations
xys = np.vstack((xx.ravel(), yy.ravel())).T
# Create a density map
x_shape = xx.shape
# Get the kernel density estimator
kde = KernelDensity(bandwidth=bandwidth,metric=metric,
kernel=kernel, algorithm='ball_tree')
# Fit it to the coordinate pairs
_ = kde.fit(positions)
# Evaluate
z = np.exp(kde.score_samples(xys))
print(np.max(z))
zi = np.arange(xys.shape[0])
# Plug densities into grid
zg = -9999 np.zeros(xys.shape[0])
zg[zi] = z
xyz = np.hstack((xys[:, :2], zg[:, None]))
# Get the density values arranged on the grid
z = xyz[:, 2].reshape(x_shape)
temp = z[::-1, :]
output_arr = temp.reshape(-1, temp.shape[0], temp.shape[1])
# Write the densities to a raster
with rasterio.open(
output_raster,
'w',
driver='GTiff',
height=output_arr.shape[1],
width=output_arr.shape[2],
dtype=output_arr.dtype,
crs=CRS.from_epsg(4326),
count=1,
transform=rasterio.transform.from_bounds(x_min, y_min, x_max, y_max, output_arr.shape[2], output_arr.shape[1])
) as dst:
dst.write(output_arr)
if __name__ == "__main__":
positions = [[126.82800884821953, 8.021550450814345],
[123.0835913004416, 15.887493017360754],
[122.87172138544588, 15.155979776107289],
[122.48465193221716, 15.233649683534475],
[122.26320643954872, 16.71625103407011],
[122.13275884500477, 15.941644592949958],
[120.63772441542471, 7.078277119741588],
[120.57180822188472, 7.537689414917545],
[119.53047809084589, 1.396741864447578],
[119.51652407635684, 1.7028166423529711],
[119.35538543402562, 7.795232293743844],
[119.35371605376332, 1.7139590065581176],
[118.21983976700818, 0.2725608428591114],
[116.32507063966972, -2.0478066628388163],
[115.9455871941716, -2.2758686356158915],
[110.54879990595637, 4.849182291868757],
[109.00373897612512, 12.330559666134512],
[108.56317006080423, 23.10356852435795],
[107.95374212609899, -3.878293744564539],
[107.6618148392204, -4.215545933851648],
[107.39598092145678, -3.3557991558597426],
[107.38347877309276, -4.243848824653475],
[107.3802332907293, -4.724984303635246],
[106.92298020128571, 3.3377440975999058],
[106.8467663232349, -3.427384435159751],
[106.6198566766759, 3.327030211530555],
[106.59035576911651, 3.409433089119516],
[106.48649132403538, 3.5936300679047966],
[106.2879431146126, 3.039670857739856],
[105.96323043582797, 2.5103916023650656],
[105.9540323861389, 2.596746532847891],
[105.80111748849575, 3.388380151516756],
[105.62119198430719, 3.2169296961449554],
[105.43276377101233, 2.6840109661437204],
[105.29236334314527, 2.420170430982717],
[104.94141265184744, 3.091707354213681],
[103.08902291491331, 3.1932135322924133],
[102.59488296531873, 14.93503092216549],
[100.7213889691745, 5.834246665586201],
[100.70491932538964, 5.2594820067014245],
[100.51665775078591, 6.0369426594605855],
[100.51156199546038, 5.491942119998682],
[100.45311457726862, 5.281343969279209],
[99.984116880471, 5.658350660638604],
[93.51170627287425, 24.024373245961645],
[93.34991893283902, 23.04050533807432],
[84.93884193888668, 19.384547030288207],
[84.30999142795147, 18.825326243832105],
[84.1630944193751, 19.06013889689632],
[83.80094785724114, 18.57306909774846],
[74.16321921976069, 23.579347585345776],
[72.4113965790803, 21.875517403679595],
[49.40472412468231, 32.2487630729451],
[42.90510332039255, -12.821849510976579],
[42.408207428324495, -12.31050970009727],
[42.36825610793828, -13.083052941231413],
[42.30285486383656, -12.234780003717532],
[15.328057669295298, -7.460883355600632],
[14.631592099379093, -7.440778982157976],
[14.563929300312948, -7.140268202440664],
[14.446656807020666, -4.699494598106393],
[14.188788859460905, -6.430418645148537],
[13.44490187975298, -2.8654279482460323],
[13.301089335672936, -2.593387816196834],
[13.131727857324034, -3.412434046655619],
[11.637624067618695, 5.306602656962694],
[11.537324701566494, 1.5773310360579327],
[11.056051828014489, 5.372994263069668],
[10.981944105212998, 6.05789466930291],
[10.978615683124655, 5.7586879077143225],
[10.384229532923067, 2.6509917300959476],
[10.293978958054748, 5.6087142487617045],
[9.724503564938162, 5.965801337392755],
[9.228154036572047, 6.4564328707855605],
[8.847083818460739, 4.696640992862242],
[8.724622829999017, 5.5476494764785516],
[8.483278678008926, 6.612624047942372],
[8.44366045716664, 6.2122982089038725],
[8.4255624128847, 4.755664077859387],
[8.11860899795907, 5.659724263701104],
[7.912362077517271, 4.87480562915889],
[7.563449250527216, 2.842579773546474],
[7.2608575851074, 5.16577516485171],
[7.004069229900638, 3.5416918941072804],
[6.9915716303567494, 3.7362296571866294],
[6.468876406999725, 5.010859767233725],
[6.203147917904825, 4.992482439632923],
[5.4017709770599325, 7.676092696459705],
[5.350100368207385, 7.762605113995827],
[5.279221956366327, 4.915935839020336],
[5.213104554080347, 8.281676925077297],
[5.1108484406102805, 7.9040681892696485],
[5.059337403465768, 8.140534352024792],
[4.861618772269268, 8.322655646328752],
[4.80376638793241, 8.062341031849334],
[4.665446704573248, 7.477404025788393],
[4.6477402888853145, 7.797020093234158],
[4.609044098910636, 38.765860093618905],
[4.555126307535386, 7.873929016757312],
[4.4195324599539845, 7.394848626095032],
[4.400283930670644, 8.038284539940614],
[4.347819621721147, 8.443859742876246],
[4.240704264765369, 6.955830447603886],
[4.227870824209585, 5.751072313355475],
[4.033821062618696, 7.0740805209122595],
[3.665972118522844, 6.545536856751896],
[3.4165849005141005, 7.191717476638518],
[3.121450235674562, 8.103710628355616],
[1.8057346437941182, 1.3314371195302515],
[0.21998421850813876, 6.744306925430884],
[-12.310298533627448, 11.362835062050264],
[-49.352317054841336, 2.010101652464972],
[-49.56587070660965, 1.366869361066606],
[-49.5821860267535, 1.824258170311353],
[-70.58665807820438, 20.03257364630837],
[-70.6803277335339, 19.902301232265422],
[-70.78620439744233, 20.024999949922996],
[-70.86459827149523, 20.273742251629713],
[-71.02033226779315, 19.891866165854587],
[-73.57317798569044, 12.265930473198331],
[-75.32300214385347, -10.734649751468147],
[-75.36631826293349, -10.206201123969526],
[-75.37463804230384, -10.724232696199014],
[-75.40829227919468, -10.817431611704407],
[-75.46984739081694, -10.195876463554633],
[-75.56266706716431, -10.202240256127965],
[-75.74233061116121, -10.647556252995775],
[-75.90503122834087, -10.297561312609464],
[-75.94114020328095, -10.530481915516726],
[-78.13302896559648, -1.2629721839381856],
[-78.42506520505198, -0.6805387090496724],
[-78.68351568134375, -1.1006283268898114],
[-79.09221180056895, -1.5423219306900116],
[-90.05839881111541, 21.022199691388156],
[-91.3208074507767, 20.58263399988673],
[-91.86906142999138, 20.169783366358622],
[-91.89838954465436, 20.49386425203851]]
bandwidth = 1.0
cell_size = 0.1
extent = [-180, 180, -90, 90]
metric = "euclidean"
kernel = "gaussian"
output_raster = metric "_" kernel "_" str(bandwidth).split(".")[0] ".tif"
# The parameters that I think should do the trick
kernel_density_lat_lon(positions, bandwidth, metric, kernel, cell_size,
extent, output_raster)
# The parameters that get me closest to the desired output
# This requires multiplying all of the density probabilities by 205...
bandwidth = 1.0
cell_size = 0.1
extent = [-180, 180, -90, 90]
metric = "euclidean"
kernel = "epanechnikov"
multiplier = 205
output_raster = metric "_" kernel "_" str(bandwidth).split(".")[0] ".tif"
kernel_density_lat_lon(positions, bandwidth, metric, kernel, cell_size,
extent, output_raster, multiplier)
Я много думал об этом, и я в тупике относительно того, почему оценки плотности меньше, чем я ожидал. Спасибо за вашу помощь.
Комментарии:
1. Интересная проблема. «преобразование моих входных данных lat / lon в радианы» <- Где это происходит?
2. Спасибо за ваш ответ. Приведенная выше реализация не преобразуется в радианы. Связанный пример делает. Когда я наткнулся на пример, я подумал, что это именно то, что мне нужно. Я заменил данные примера своими собственными, но получил забавные результаты, в частности, над центром 2d-решетки появился большой эллипс ячеек с меньшим значением. Приведенная выше реализация использует мое текущее лучшее мышление: расположите точки широты / широты на декартовой плоскости и используйте евклидову метрику. Реальная проблема заключается в том, почему значения плотности такие низкие.
3. Широта и долгота — это сферическая система координат, описывающая местоположения на поверхности 3D-сферы, а не 2D-плоскости. Таким образом, ваша плотность ядра 2D здесь не имеет смысла. Координаты должны быть в картографической проекции (декартова система координат 2D).
4. Вы также пытались преобразовать свою пропускную способность в радианы? Поскольку широта / долгота указаны в радианах, пропускная способность также должна быть в радианах.