Введение
Начну с конца. Это скриншот с некой web-карты, визуализирующей среднюю стоимость недвижимости на вторичном рынке Саратова и Энгельса:
Цвета на карте можно соотнести с цветами на «легенде», цвет на «легенде» соответствует средней стоимости квадратного метра общей площади в тысячах рублей.
Точка на карте соответствует одному предложению по продаже (на вторичном рынке) квартиры с Авито. Всего таких точек, как видно на «легенде», для построения графика использовалось 4943.
Карта в интерактивном виде доступна на GitHub.
А теперь немного предыстории..
Давным-давно…
… когда будучи студентом мехмата я писал работу по статистической обработке данных по ценам на недвижимость, я случайно нашел полезный документ, — градация районов по ценам на недвижимость, что-то вроде “Ценовые зоны Саратова”. Не знаю, кто ее составлял, но составлена она была правдоподобно, — анализ показал, что (если учесть другие факторы), разница между средней стоимостью квадратного метра квартир из соседних ценовых зон составляла примерно 10%. С тех выступать в роли оценщика недвижимости (на бытовом уровне), как и большинству взрослых людей, приходилось не раз.
Идея снова подойти к этому вопросу более “научно” пришла в голову когда я увидел публикацию Jeff Kaufman Boston Apartment Price Maps.
Товарищ сделал следующее:
- Python скрипт, который парсит с сайта PadMapper данные по стоимости аренды недвижимости.
- еще скрипт, который используя методы Spatial interpolation, визуализирует эти данные в виде растрового изображения.
- Наложил этот рисунок на страничку с Google Maps
(Исходники и данные — здесь)
В результате у него получилась такая интерактивная карта распределения стоимости аренды на карте Бостон-а:
После первого “Вау, прикольно”, следующая мысль — “А что если попробовать сделать такое же, но для России, только не по стоимости аренды, а по ценам на квадрат недвижимости”.
Сразу оговорюсь, я, по текущему роду деятельности, скорее database developer, и не являюсь специалистом ни по GIS данным, ни по front-end разработке, ни по статистике, ни по Python, поэтому отнеситесь снисходительно к… еще не знаю к чему, но надеюсь этого будет не очень много и вы мне об этом напишете для исправления.
И здесь возникает первая крупная задача
Web Scraping с Авито.
Так как с Python-ом на тот момент я был совсем на “вы”, после изучения HTML кода страниц Авито, свой собственный велосипед я решил написать в виде библиотеки на C#, используя для парсинга AngleSharp.
Не буду сейчас углубляться в особенности реализации, замечу только что:
- все метаданные парсинга Авито я вынес в отдельный конфигурационный класс, чтобы в будущем, при изменении структуры разметки на Авито, можно было бы просто изменить строку в конфигурационном файле.
- дабы избежать бана от Авито, пришлось добавлять искусственные задержки в 10 секунд после обработки каждого объявления. (Согласен, — не самый умный способ, но найти подходящий пул прокси -айпишников навскидку не удалось).
В конечном итоге, запустив парсилку на ночь, я получил нужные данные, записанные в табличку в базе
Очистив данные от “Старого фонда”, “Элитных”, и всяких явно неадекватных предложений, я получил выборку из пяти тысяч точек — продажи квартир на вторичном рынке Саратова и Энгельса. И это дало возможность перейти к следующему шагу:
Обработка GIS данных.
Общая постановка задачи:
— есть дискретный набор значений некой метрики(температура, высота, концентрация вещества и т.п), с привязкой к координатам на плоскости (X, Y или широта, долгота),
— нужно иметь возможность предсказать значение этой метрики в произвольной точке этой плоскости,
— результат визуализируется в виде цветовой схемы на плоскости. Для цветовой схемы для удобства можно дискретизировать цвета.
Пример, — сняли в один момент показания с датчиков температуры в разных местах, обработали, получили такую картинку:
На самом деле, задачи такого рода требуют погружения в тему “Spatial interpolation”. Здесь есть свой математический аппарат, например:
— Inverse Distance Weighting (IDW) Interpolation
— Kriging
И есть свой инструментарий в составе популярных GIS пакетов, например
— Интерполяция точечных данных qGIS
— GeoStatistical Analyst ArcGIS — — Use of SAGA GIS for spatial interpolation(kriging)
— Kriging in R.
IDW в составе qGIS (см первую из ссылок) я попробовал, — результат меня не удовлетворил.
Поэтому, чтобы не завязнуть надолго в освоении нового инструментария, я предпочел воспользоваться готовым Рython скриптом от Jeff Kraufman . Для расчета расстояний и для трансляции GPS координат в X,Y координаты на растре здесь применяется упрощенная линейная формула, без применения используемой в Google MapsWeb Mercator проекции.
Насколько я понял, в скрипте используется вариация на тему 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 карты Бостона от Jeff Kraufman. Немного видоизменяем его под свои данные, и получаем вот такую «картину маслом».
Для проверки на явные “косяки” в Georeferencing загружаем тот же набор точек в qGIS, подгружаем ГуглоКарты, сравниваем картинки.
Убедившись, что точки на нарисованной карте находятся на своих местах, а цвета на ней правдоподобно отражают цены на квартиры в нашей деревне, можем выдохнуть — “Yes! Работает!”.
Только вот немного огорчает, что картинка у нас залезла в Волгу. Что естественно, — алгоритм расчета ничего про Волгу не знает.
Стало быть возникает задача:
Обрезать изображение на карте по береговой линии
Попытка найти гео-данные по административным границам городов с учетом береговой линии не увенчалась успехом. Есть Саратов, есть Энгельс, и граница между ними проходит посреди Волги. Ладно, мы пойдем другим путем.
Скачиваем из открытых источников shape file c данными по береговым линиям, загружаем как векторный слой в qGIS, выбираем нужные полигоны, разбираемся в структуре файлов, вытаскиваем полигоны, которые относятся к Волге у Саратова, заливаем данные в БД.
Немного магии типа:
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 у меня отношения не сложились, и после нескольких попыток я решил, что используя 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.
Получаем API_KEY, выясняем, что:
— максимальное разрешение для бесплатного использования =640*640
— параметризации по BoundingBox -, в отличие от BING static maps API гугл не поддерживает, и придется дополнительно соотносить два рисунка в координатной плоскости.
Ладно, пробуем. Получаем от 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 в папочку, расшаренную для , и можем хвастаться перед IT-сообществом.
И что дальше?
А дальше возникает много вопросов, на которые я пока не знаю четкого ответа:
- Насколько корректно работает алгоритм расчета здесь по сравнению с другими методами spatial interpolations
- Какие еще факторы, кроме расположения объекта недвижимости, имеет смысл учитывать при расчете и как ввести в расчетную модель
- Можно ли сделать постоянное обновление данных в базе и построить на этом какую то систему оценку недвижимости. (Если меня читают на этом сайте представители Авито, хотелось бы услышать комментарии)
Ну и… спасибо всем дочитавшим до этого места. Мне было интересно решить практическую задачу, которая оказалась, “rather challenging for me”, надеюсь вам было интересно прочитать о моем опыте.
Автор: dklmn
