- PVSM.RU - https://www.pvsm.ru -
Начну с конца. Это скриншот с некой web-карты, визуализирующей среднюю стоимость недвижимости на вторичном рынке Саратова и Энгельса:
Цвета на карте можно соотнести с цветами на «легенде», цвет на «легенде» соответствует средней стоимости квадратного метра общей площади в тысячах рублей.
Точка на карте соответствует одному предложению по продаже (на вторичном рынке) квартиры с Авито. Всего таких точек, как видно на «легенде», для построения графика использовалось 4943.
Карта в интерактивном виде доступна на GitHub [1].
Давным-давно…
… когда будучи студентом мехмата я писал работу по статистической обработке данных по ценам на недвижимость, я случайно нашел полезный документ, — градация районов по ценам на недвижимость, что-то вроде “Ценовые зоны Саратова”. Не знаю, кто ее составлял, но составлена она была правдоподобно, — анализ показал, что (если учесть другие факторы), разница между средней стоимостью квадратного метра квартир из соседних ценовых зон составляла примерно 10%. С тех выступать в роли оценщика недвижимости (на бытовом уровне), как и большинству взрослых людей, приходилось не раз.
Идея снова подойти к этому вопросу более “научно” пришла в голову когда я увидел публикацию Jeff Kaufman Boston Apartment Price Maps [2].
Товарищ сделал следующее:
(Исходники и данные — здесь [5])
В результате у него получилась такая интерактивная карта распределения стоимости аренды на карте Бостон [6]-а:
После первого “Вау, прикольно”, следующая мысль — “А что если попробовать сделать такое же, но для России, только не по стоимости аренды, а по ценам на квадрат недвижимости”.
Сразу оговорюсь, я, по текущему роду деятельности, скорее database developer, и не являюсь специалистом ни по GIS данным, ни по front-end разработке, ни по статистике, ни по Python, поэтому отнеситесь снисходительно к… еще не знаю к чему, но надеюсь этого будет не очень много и вы мне об этом напишете для исправления.
И здесь возникает первая крупная задача
Так как с Python-ом на тот момент я был совсем на “вы”, после изучения HTML кода страниц Авито, свой собственный велосипед я решил написать в виде библиотеки на C#, используя для парсинга AngleSharp [7].
Не буду сейчас углубляться в особенности реализации, замечу только что:
В конечном итоге, запустив парсилку на ночь, я получил нужные данные, записанные в табличку в базе
Очистив данные от “Старого фонда”, “Элитных”, и всяких явно неадекватных предложений, я получил выборку из пяти тысяч точек — продажи квартир на вторичном рынке Саратова и Энгельса. И это дало возможность перейти к следующему шагу:
Общая постановка задачи:
— есть дискретный набор значений некой метрики(температура, высота, концентрация вещества и т.п), с привязкой к координатам на плоскости (X, Y или широта, долгота),
— нужно иметь возможность предсказать значение этой метрики в произвольной точке этой плоскости,
— результат визуализируется в виде цветовой схемы на плоскости. Для цветовой схемы для удобства можно дискретизировать цвета.
Пример, — сняли в один момент показания с датчиков температуры в разных местах, обработали, получили такую картинку:
На самом деле, задачи такого рода требуют погружения в тему “Spatial interpolation [8]”. Здесь есть свой математический аппарат, например:
— Inverse Distance Weighting (IDW) Interpolation [9]
— Kriging [10]
И есть свой инструментарий в составе популярных GIS пакетов, например
— Интерполяция точечных данных qGIS [11]
— GeoStatistical Analyst ArcGIS [12] — — Use of SAGA GIS for spatial interpolation(kriging) [13]
— Kriging in R [14].
IDW в составе qGIS (см первую из ссылок) я попробовал, — результат меня не удовлетворил.
Поэтому, чтобы не завязнуть надолго в освоении нового инструментария, я предпочел воспользоваться готовым Рython скриптом от Jeff Kraufman [15]. Для расчета расстояний и для трансляции GPS координат в X,Y координаты на растре здесь применяется упрощенная линейная формула, без применения используемой в Google MapsWeb Mercator проекции [16].
Насколько я понял, в скрипте используется вариация на тему Inverse Distance Weighting (IDW), главная (и самая тяжеловесная по расчетной части) в этом скрипте часть здесь в этом коде.
gaussian_variance = IGNORE_DIST/2
gaussian_a = 1 / (gaussian_variance * math.sqrt(2 * math.pi))
gaussian_negative_inverse_twice_variance_squared = -1 / (2 * gaussian_variance * gaussian_variance)
def gaussian(prices, lat, lon, ignore=None):
num = 0
dnm = 0
c = 0
for price, plat, plon, _ in prices:
if ignore:
ilat, ilon = ignore
if distance_squared(plat, plon, ilat, ilon) < 0.0001:
continue
weight = gaussian_a * math.exp(distance_squared(lat,lon,plat,plon) *
gaussian_negative_inverse_twice_variance_squared)
num += price * weight
dnm += weight
if weight > 2:
c += 1
# don't display any averages that don't take into account at least five data points with significant weight
if c < 5:
return None
return num/dnm
def distance_squared(x1,y1,x2,y2):
return (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)
Подогнал свой массив данных под нужный формат, убрал из скрипта ту часть, которая делает расчет корреляций относительно bedroom_category (количество комнат), поигрался с входными параметрами, запустил расчет. На моем компе (Intel® Core(TM) i7-6700 CPU, 16GB RAM), расчет с 5 тыс точек на входе и разрешением изображения на выходе 1000*1000 занял примерно полчаса.
Осталось немного поработать с Front-End -ом. Смотрим HTML source карты Бостона [17]от Jeff Kraufman. Немного видоизменяем его под свои данные, и получаем вот такую «картину маслом».
Для проверки на явные “косяки” в Georeferencing [18] загружаем тот же набор точек в qGIS, подгружаем ГуглоКарты, сравниваем картинки.
Убедившись, что точки на нарисованной карте находятся на своих местах, а цвета на ней правдоподобно отражают цены на квартиры в нашей деревне, можем выдохнуть — “Yes! Работает!”.
Только вот немного огорчает, что картинка у нас залезла в Волгу. Что естественно, — алгоритм расчета ничего про Волгу не знает.
Стало быть возникает задача:
Попытка найти гео-данные по административным границам городов с учетом береговой линии не увенчалась успехом. Есть Саратов, есть Энгельс, и граница между ними проходит посреди Волги. Ладно, мы пойдем другим путем.
Скачиваем из открытых источников shape [19]file c данными по береговым линиям, загружаем как векторный слой в qGIS, выбираем нужные полигоны, разбираемся в структуре файлов, вытаскиваем полигоны, которые относятся к Волге у Саратова, заливаем данные в БД.
Немного магии [20] типа:
DECLARE @Volga geography;
select @Volga= PlacePolygon.MakeValid() from PlacesPolygons where PlacesPolygonID =1003
DECLARE @Saratov geography;
select @Saratov= PlacePolygon.MakeValid() from PlacesPolygons where PlacesPolygonID =2
select @Saratov.STDifference(@Volga)
И мы получаем нужный полигон — границы города с учетом береговой линии.
Остается загрузить это в qGIS, загрузить туда же сгенерированную скриптом картинку как растровый слой, провести Georeferencing для этой картинки, и обрезать ее по полигону с границами городов, все стандартными средства qGIS. Все было бы здорово если бы… мне это удалось сделать. Но увы… с Georeferencing в qGIS [21] у меня отношения не сложились, и после нескольких попыток я решил, что используя spatial функции из Microsoft.SqlServer.Types, проще быстренько сделать очередной велосипед в виде утилиты на шарпе.
string GPSPolygonVolga = " POLYGON ((46.2324447632 51…. ))";
private void CropByGPS()
{
Bitmap bmp = new Bitmap(inputFilePath);
int w = bmp.Width;
int h = bmp.Height;
var GpsSizeOfPyxelX = (MAX_LON - MIN_LON) / w;
var GpsSizeOfPyxelY = (MAX_LAT - MIN_LAT) / h;
var VolgaPolygon = SqlGeography.STPolyFromText(new System.Data.SqlTypes.SqlChars(GPSPolygonVolga), SRID);
for (int y = 0; y < h; y++)
{
for (int x = 0; x < w; x++)
{
Color c = bmp.GetPixel(x, y);
var GPS = SqlGeography.Point(MAX_LAT - GpsSizeOfPyxelY * y, MIN_LON + GpsSizeOfPyxelX * x, SRID);
if (VolgaPolygon.STContains(GPS)) //попадает ли точка в "водяной" полигон
bmp.SetPixel(x, y, Color.Transparent);
}
}
bmp.Save(outputFilePath, System.Drawing.Imaging.ImageFormat.Png);
}
Загрузили исходный файл с растром, получили на выходе тот же файл, но в немного “обкусанном виде”. Наложили на Google Maps в браузере, смотрим, — на Волгу ничего не залезает, отлично! Смотрим чуть внимательнее и видим… что кажется кое-где откусили лишнее. Косяк в данных или косяк в процедуре? Снова смотрим на загруженные слои в qGIS.
Да… действительно, кое-где у нас порой… полигоны из слоя водоема явно залезает в город. Прям не Энгельс, а Венеция какая-то… Ладно, нет стремления к совершенству, мы пойдем другим путем.
Раз векторные данные вышли из доверия, будем работать с растровыми. Нам нужно наложить две картинки, и вырезать из первой те точки, которые соответствуют “голубому” во второй.
Картинку с картой нам могут дать, например Google Maps static API [22].
Получаем API_KEY, выясняем, что:
— максимальное разрешение для бесплатного использования =640*640
— параметризации по BoundingBox [23]-, в отличие от BING static maps API [24] гугл не поддерживает, и придется дополнительно соотносить два рисунка в координатной плоскости.
Ладно, пробуем. Получаем от Google Statics Map API некую картинку с картой, покрывающую интересующую нас область, пишем такой код:
private static void CropByMapImage()
{
Bitmap bmp = new Bitmap(inputFilePath);
Bitmap bmpMap = new Bitmap(mapToCompareFilePath);
int w = bmp.Width;
int h = bmp.Height;
var GpsSizeOfPyxelX = (MAX_LON - MIN_LON) / w;
var GpsSizeOfPyxelY = (MAX_LAT - MIN_LAT) / h;
int wm = bmpMap.Width;
int hm = bmpMap.Height;
var cntAll = w * h;
var cntRiver = 0;
for (int y = 0; y < h; y++)
{
for (int x = 0; x < w; x++)
{
var currentLatitude = MAX_LAT - GpsSizeOfPyxelY * y;
var currentLongitude = MIN_LON + GpsSizeOfPyxelX * x;
if (currentLongitude < MIN_LON_M || currentLongitude > MAX_LON_M || currentLatitude < MIN_LAT_M || currentLatitude > MAX_LAT_M) continue; //todo - something with this ugly code
var CorrespondentX = (int)((currentLongitude - MIN_LON_M) / (MAX_LON_M - MIN_LON_M) * wm);
var CorrespondentY = (int)((MAX_LAT_M-currentLatitude) / (MAX_LAT_M - MIN_LAT_M) * hm);
Color cp = bmp.GetPixel(x, y); //color of current pixel
Color mp = bmpMap.GetPixel(CorrespondentX, CorrespondentY); //color of correspondent Map pixel
if (Math.Max(Math.Max( Math.Abs(mp.R - WaterMapColor.R) , Math.Abs(mp.G - WaterMapColor.G)) , Math.Abs(mp.B -WaterMapColor.B))< colorDiff /* && cp.A != Color.Transparent.A */)// blue water is not always the same blue
{
bmp.SetPixel(x, y, Color.Transparent);
cntRiver++;
}
}
}
bmp.Save(outputFilePath, System.Drawing.Imaging.ImageFormat.Png);
}
Обрабатываем наш файлик с разрешением 1000*1000 с помощью файлика с разрешением 640*640, получаем… такой результат (здесь для наглядности те точки, которые должны быть Transparent, сделаны желтыми):
Результат — печальный, в том числе еще и потому, что Гуглокарта масштабируется дискретно, и поэтому с большим запасом перекрывает наш сгенерированный растр. То есть или нужно брать для “обрезки” много карт с подходящим масштабом или… как то все таки найти одну, но с хорошим разрешением.
Я — пошел вторым путем, то есть просто за-скриншотил нужную проекцию обычной гуглокарты из браузера. GPS координаты углов (mapping bound), я узнал, добавив в JavaScript code такое:
google.maps.event.addListener(map, 'bounds_changed', function () {
console.log(map.getBounds());
Обработал наш растр с помощью новой карты из сохраненного скриншота, загрузил в браузер и… результат меня — скорее порадовал. Достаточно корректно(в пределах точности GPS позиционирования) “зачистился” даже маленький пруд в местном городском парке.
Все, будем считать первый этап квеста завершенным, заливаем результаты на GitHub в папочку, расшаренную для
А дальше возникает много вопросов, на которые я пока не знаю четкого ответа:
Ну и… спасибо всем дочитавшим до этого места. Мне было интересно решить практическую задачу, которая оказалась, “rather challenging for me”, надеюсь вам было интересно прочитать о моем опыте.
Автор: dklmn
Источник [26]
Сайт-источник PVSM.RU: https://www.pvsm.ru
Путь до страницы источника: https://www.pvsm.ru/c-2/250460
Ссылки в тексте:
[1] GitHub: https://dmitriik.github.io/RealtyEstimation/
[2] Boston Apartment Price Maps: https://www.jefftk.com/p/updated-boston-apartment-price-maps
[3] PadMapper: https://www.padmapper.com
[4] на страничку с Google Maps: https://www.jefftk.com/apartment_prices/
[5] здесь: https://github.com/jeffkaufman/apartment_prices
[6] Бостон: https://www.youtube.com/watch?v=x-64CaD8GXw
[7] AngleSharp: https://anglesharp.github.io/
[8] Spatial interpolation: http://www.bisolutions.us/A-Brief-Introduction-to-Spatial-Interpolation.php
[9] Inverse Distance Weighting (IDW) Interpolation : http://http//gisgeography.com/inverse-distance-weighting-idw-interpolation/
[10] Kriging: https://en.wikipedia.org/wiki/Kriging
[11] Интерполяция точечных данных qGIS: http://www.qgistutorials.com/ru/docs/interpolating_point_data.html
[12] GeoStatistical Analyst ArcGIS : http://npk-kaluga.ru/GeoStatAnalyze_AGIS.htm
[13] Use of SAGA GIS for spatial interpolation(kriging) : http://www.dmcsee.org/uploads/file/337_4_saga_kriging_manual.pdf
[14] Kriging in R: https://rpubs.com/nabilabd/118172
[15] готовым Рython скриптом от Jeff Kraufman : http://https//github.com/jeffkaufman/apartment_prices
[16] используемой в Google MapsWeb Mercator проекции: https://en.wikipedia.org/wiki/Web_Mercator
[17] HTML source карты Бостона : http://HTML source карты Бостона
[18] Georeferencing: https://en.wikipedia.org/wiki/Georeferencing
[19] shape : https://ru.wikipedia.org/wiki/Shapefile
[20] магии: https://msdn.microsoft.com/en-us/library/cc280766.aspx
[21] Georeferencing в qGIS: http://www.qgistutorials.com/en/docs/advanced_georeferencing.html
[22] Google Maps static API: http://https//developers.google.com/maps/documentation/static-maps/?hl=ru
[23] BoundingBox : http://wiki.openstreetmap.org/wiki/Bounding_Box
[24] BING static maps API: https://msdn.microsoft.com/en-us/library/ff701724.aspx
[25] хостинга: https://www.reg.ru/?rlink=reflink-717
[26] Источник: https://habrahabr.ru/post/324596/?utm_source=habrahabr&utm_medium=rss&utm_campaign=sandbox
Нажмите здесь для печати.