Тема
Правила форума
Для предотврашения спама первые сообщения вновь зарегистрированных пользователей проходят ручную премодерацию.
Для предотврашения спама первые сообщения вновь зарегистрированных пользователей проходят ручную премодерацию.
Тема
Киньте ссылкой на теорию, плизз :)
Есть N точек с растровыми координатами. Это N ооочень велико :)
Знаю датум, знаю проекцию, знаю геодезические координаты нескольки таких точек.
Нужно найти такие координаты для всех остальных. Тоесть сделать то, что умеет Ози после привязки растра. Насколько я понял, нужно по известным данным посчитать параметры аффинного преобразования, а потом его выполнить. И что-то еще...
Есть N точек с растровыми координатами. Это N ооочень велико :)
Знаю датум, знаю проекцию, знаю геодезические координаты нескольки таких точек.
Нужно найти такие координаты для всех остальных. Тоесть сделать то, что умеет Ози после привязки растра. Насколько я понял, нужно по известным данным посчитать параметры аффинного преобразования, а потом его выполнить. И что-то еще...
|
||
-
- Сообщения: 18
- Зарегистрирован: 27 фев 2007, 21:17
Re: Тема
Читай тут http://www-ec.njit.edu/surveying/Docume ... .pdf но если задача однократная, то проще сделать это какой-нибудь ГИС, это многие умеют (рекомендую GeoMedia, ArcGIS)
Re: Геопривязка своими руками
Задача однократная, но растра как такового нету, есть просто куча координат в базе, которые можно считать растровыми. Конечно, можно по ним сформировать картинку. Но, думаю, удобнее будет вывести формулы преобразования и прогнать через них все мои точки, тем более, что их сразу же надо положить в другую базу.
Re: Тема
Спасибо за ссылку. Полезная информация
Re: Геопривязка своими руками
Если умеете немного программировать, то задача довольно простая.
Пишете линейные уравнения (или квадратичные, если есть искажения), связывающие координаты растра и топо (типа: X=a*x+b*y+m; Y=c*x+d*y+n). Будете иметь 6 неизвестных параметров (a,b,m,c,d,n). Ошибка - сумма квадратов расхождений по X и Y для всех точек с известными координатами. Берёте частные производные от ошибки по всем неизвестным параметрам и приравниваете 0 (ищете минимум квадратичной ошибки)- получите 2 системы линейных уравнений, каждая с 3 неизвестными. Решение - неизвестные параметры. Прогоняете файл со своими точками через эти линейные уравнения и пишете в файл. Понятно, при нелинейной интерполяции параметров будет больше.
Пишете линейные уравнения (или квадратичные, если есть искажения), связывающие координаты растра и топо (типа: X=a*x+b*y+m; Y=c*x+d*y+n). Будете иметь 6 неизвестных параметров (a,b,m,c,d,n). Ошибка - сумма квадратов расхождений по X и Y для всех точек с известными координатами. Берёте частные производные от ошибки по всем неизвестным параметрам и приравниваете 0 (ищете минимум квадратичной ошибки)- получите 2 системы линейных уравнений, каждая с 3 неизвестными. Решение - неизвестные параметры. Прогоняете файл со своими точками через эти линейные уравнения и пишете в файл. Понятно, при нелинейной интерполяции параметров будет больше.
|
||
Re: Геопривязка своими руками
С программированием не расстаюсь уже лет 6 :). Вспомнилась лаба по интерполяции, что мы делали года 3 назад. Кстати, для преобразования растровых координат в координаты проекции Меркатора на сферу нужна нелинейная интерполяция?
Re: Геопривязка своими руками
Если точки покрывают континенты, то лучше ждать ответа от монстров форума. На маленьких площадях вся меркаторность и эллипсоидальность ИМХО уже учтена в координатах известных точек. Впрочем, даже для континентов, если точки с известными координатами более или менее равномерно покрывают интересующую площадь, можно разбить всю область на небольшие куски и выполнить линейное преобразование внутри каждого из них. Т.е. прокластеризовать пространство вокруг известных точек.
Re: Геопривязка своими руками
: Есть N точек с растровыми координатами. Это N ооочень велико :)
: Знаю датум, знаю проекцию, знаю геодезические координаты нескольки таких точек.
: Нужно найти такие координаты для всех остальных. Тоесть сделать то, что умеет Ози после привязки растра. Насколько я понял, нужно по известным данным посчитать параметры аффинного преобразования, а потом его выполнить. И что-то еще...
1. Пересчитать географические координаты точек привязки в датум и координаты проекции (которые известны).
2. Рассчитать параметры афинного преобразования по известным парам координат ("растровых" и координат проекции) для точек привязки.
3. При помощи афинного преобразования рассчитать координаты проекции для остальных точек.
4. Пересчитать координаты проеции для остальных точек в географические координаты и соответствующий датум.
: Знаю датум, знаю проекцию, знаю геодезические координаты нескольки таких точек.
: Нужно найти такие координаты для всех остальных. Тоесть сделать то, что умеет Ози после привязки растра. Насколько я понял, нужно по известным данным посчитать параметры аффинного преобразования, а потом его выполнить. И что-то еще...
1. Пересчитать географические координаты точек привязки в датум и координаты проекции (которые известны).
2. Рассчитать параметры афинного преобразования по известным парам координат ("растровых" и координат проекции) для точек привязки.
3. При помощи афинного преобразования рассчитать координаты проекции для остальных точек.
4. Пересчитать координаты проеции для остальных точек в географические координаты и соответствующий датум.
Не знаю может поможет..
Вот надысь сделал аддон к мапедиту где преобразовывал координаты из WGS в UTM и обратно..
Ниже привожу исходник оной функции на C++
Эта функция преобразует координаты в угловых координатах(в градусах) в координаты проекции UTM и возвращает значение Центрального меридиана зоны (в градусах).
int WGStoUTM (float *Lat, float *Lon)
{
int CentrMeridian;
float FE=500000; //false easting
float FN=0; //false northern
CentrMeridian=GetCM(*Lon); //central meridian of zone
float a=6378137; //semi-major axis of ellipsoid in WGS84
float f=1/298.257223; //flattening
float K0=0.9996013; //distorsion coefficient
float LatRad=*Lat/ 57.29577951; //point latitude in radians
float LonRad=*Lon/ 57.29577951; //point longitude in radians
float Lat0=0/57.29577951; //latitude of origin
float Lon0=CentrMeridian/ 57.29577951; //longitude of origin in radians
float esq=(2*f-pow(f,2));
float e=sqrt(esq); //eccentricity
float e_strih=sqrt(pow(e,2)/(1-pow(e,2))); //second eccentricity
float C=esq * pow(cos(LatRad),2) / (1 - esq);
float nu=a / sqrt(1 - esq * pow(sin(LatRad),2));
float A_ = (LonRad - Lon0) * cos(LatRad);
float T = pow(tan(LatRad),2);
float E4 =pow(e,4), E6 =pow(e,6);
float M = a*(
(1 - esq/4 - 3*E4/64 - 5*E6/256 )*LatRad-
( 3*esq/8 + 3*E4/32 + 45*E6/1024)*sin(2*LatRad)+
( 15*E4/256 + 45*E6/1024 )*sin(4*LatRad)-
( 35*E6/3072 )*sin(6*LatRad)
);
float M0 = a*(
(1 - esq/4 - 3*E4/64 - 5*E6/256 )*Lat0-
( 3*esq/8 + 3*E4/32 + 45*E6/1024)*sin(2*Lat0)+
( 15*E4/256 + 45*E6/1024 )*sin(4*Lat0)-
( 35*E6/3072 )*sin(6*Lat0)
);
*Lon= FE + K0 * nu * (A_ + (1 - T + C) * pow(A_,3) / 6 + (5 - 18 * T + pow(T,2) + 72 * C - 58 * e_strih) * pow(A_,5) / 120);
*Lat= FN + K0 * (M - M0 + nu * tan(LatRad) * (pow(A_,2) / 2 + (5 - T + 9 * C + 4 * pow(C,2)) * pow(A_,4) / 24 + (61 - 58 * T + pow(T,2) + 600 * C - 330 * e_strih) * pow(A_,6) / 720));
return CentrMeridian;
}
//----------------------------------------------------------------------------------
Следующая функция принимает UTM координаты и значение Центрального меридиана зоны и преобразует их к координатам WGS84 (в градусах). Здесь прошу обратить внимание на небольшой мой косячек, в конце вычислений из значения координаты широты вычитается 0.00037. Дело в том что при обратном преобразовании координаты почему-то уходили на север на это значение. Пришлось его отнять.
Видимо где-то я что-то не так сделал, но разбираться времени небыло, оставил так.
void UTMtoWGS(float *Nord, float *East, int CentrMer)
{
float FE=500000; //false easting
float FN=0; //false northern
float a=6378137; //semi-major axis of ellipsoid in WGS84
float f=1/298.257223; //flattening
float K0=0.9996013; //distorsion coefficient
float esq=(2*f-pow(f,2)); //eccentricity in square
float e=sqrt(esq); //eccentricity
float e_strih=sqrt(pow(e,2)/(1-pow(e,2))); //second eccentricity
float e_strih_sq=pow(e,2)/(1-pow(e,2)); //second eccentricity in square
float e1=(1-pow(1-esq,0.5))/(1+pow(1-esq,0.5));
float Lat0=0/57.29577951; //latitude of origin
float E4 =pow(e,4), E6 =pow(e,6);
float Lon0=CentrMer/57.29577951;
float M0 = a*((1 - esq/4 - 3*E4/64 - 5*E6/256 )*Lat0-( 3*esq/8 + 3*E4/32 + 45*E6/1024)*sin(2*Lat0)+
( 15*E4/256 + 45*E6/1024 )*sin(4*Lat0)-( 35*E6/3072 )*sin(6*Lat0)
);
float M1=M0+(*Nord-FN)/K0;
float mu1=M1/(a*(1-esq/4-3*E4/64-5*E6/256));
float Lat1=mu1+(3*e1/2 - 27*pow(e1,3)/32)*sin(2*mu1)+
(21*pow(e1,2)/16 -55*pow(e1,4)/32)*sin(2*mu1)+
(151*pow(e1,3)/96)*sin(6*mu1)+(1097*pow(e1,4)/512)*sin(8*mu1);
float T1=tan(Lat1);
float C1=e_strih_sq*pow(cos(Lat1),2);
float nu1=a / sqrt(1 - esq * pow(sin(Lat1),2));
float D=(*East- FE)/(nu1*K0);
float ro1=a*(1-esq)/pow((1-esq*pow(sin(Lat1),2)),1.5);
*Nord=Lat1-(nu1*tan(Lat1)/ro1)*(pow(D,2)/2-(5+3*T1+10*C1-4*pow(C1,2)-9*e_strih_sq)*pow(D,4)/24+
(61+90*T1+298*C1+45*pow(T1,2)-252*e_strih_sq-3*pow(C1,2))*pow(D,6)/720 );
*East=Lon0+(D-(1+2*T1+C1)*pow(D,3)/6 + (5-2*C1+28*T1-3*pow(C1,2)+8*e_strih_sq+24*pow(T1,2))*pow(D,5)/120)/
cos(Lat1);
*Nord=*Nord*57.29577951-0.00037;
*East=*East*57.29577951;
}
//------------------------------------------------------------------------------
Ну и самое простое. Функция которая выдает значение Центрального меридиана зоны.
Ей даешь значение долготы, в градусах, она возвращает значение ЦМ.
int GetCM(size_t longitude)
{
int Lon0=(longitude/6)*6+3;
return Lon0;
}
Ниже привожу исходник оной функции на C++
Эта функция преобразует координаты в угловых координатах(в градусах) в координаты проекции UTM и возвращает значение Центрального меридиана зоны (в градусах).
int WGStoUTM (float *Lat, float *Lon)
{
int CentrMeridian;
float FE=500000; //false easting
float FN=0; //false northern
CentrMeridian=GetCM(*Lon); //central meridian of zone
float a=6378137; //semi-major axis of ellipsoid in WGS84
float f=1/298.257223; //flattening
float K0=0.9996013; //distorsion coefficient
float LatRad=*Lat/ 57.29577951; //point latitude in radians
float LonRad=*Lon/ 57.29577951; //point longitude in radians
float Lat0=0/57.29577951; //latitude of origin
float Lon0=CentrMeridian/ 57.29577951; //longitude of origin in radians
float esq=(2*f-pow(f,2));
float e=sqrt(esq); //eccentricity
float e_strih=sqrt(pow(e,2)/(1-pow(e,2))); //second eccentricity
float C=esq * pow(cos(LatRad),2) / (1 - esq);
float nu=a / sqrt(1 - esq * pow(sin(LatRad),2));
float A_ = (LonRad - Lon0) * cos(LatRad);
float T = pow(tan(LatRad),2);
float E4 =pow(e,4), E6 =pow(e,6);
float M = a*(
(1 - esq/4 - 3*E4/64 - 5*E6/256 )*LatRad-
( 3*esq/8 + 3*E4/32 + 45*E6/1024)*sin(2*LatRad)+
( 15*E4/256 + 45*E6/1024 )*sin(4*LatRad)-
( 35*E6/3072 )*sin(6*LatRad)
);
float M0 = a*(
(1 - esq/4 - 3*E4/64 - 5*E6/256 )*Lat0-
( 3*esq/8 + 3*E4/32 + 45*E6/1024)*sin(2*Lat0)+
( 15*E4/256 + 45*E6/1024 )*sin(4*Lat0)-
( 35*E6/3072 )*sin(6*Lat0)
);
*Lon= FE + K0 * nu * (A_ + (1 - T + C) * pow(A_,3) / 6 + (5 - 18 * T + pow(T,2) + 72 * C - 58 * e_strih) * pow(A_,5) / 120);
*Lat= FN + K0 * (M - M0 + nu * tan(LatRad) * (pow(A_,2) / 2 + (5 - T + 9 * C + 4 * pow(C,2)) * pow(A_,4) / 24 + (61 - 58 * T + pow(T,2) + 600 * C - 330 * e_strih) * pow(A_,6) / 720));
return CentrMeridian;
}
//----------------------------------------------------------------------------------
Следующая функция принимает UTM координаты и значение Центрального меридиана зоны и преобразует их к координатам WGS84 (в градусах). Здесь прошу обратить внимание на небольшой мой косячек, в конце вычислений из значения координаты широты вычитается 0.00037. Дело в том что при обратном преобразовании координаты почему-то уходили на север на это значение. Пришлось его отнять.
Видимо где-то я что-то не так сделал, но разбираться времени небыло, оставил так.
void UTMtoWGS(float *Nord, float *East, int CentrMer)
{
float FE=500000; //false easting
float FN=0; //false northern
float a=6378137; //semi-major axis of ellipsoid in WGS84
float f=1/298.257223; //flattening
float K0=0.9996013; //distorsion coefficient
float esq=(2*f-pow(f,2)); //eccentricity in square
float e=sqrt(esq); //eccentricity
float e_strih=sqrt(pow(e,2)/(1-pow(e,2))); //second eccentricity
float e_strih_sq=pow(e,2)/(1-pow(e,2)); //second eccentricity in square
float e1=(1-pow(1-esq,0.5))/(1+pow(1-esq,0.5));
float Lat0=0/57.29577951; //latitude of origin
float E4 =pow(e,4), E6 =pow(e,6);
float Lon0=CentrMer/57.29577951;
float M0 = a*((1 - esq/4 - 3*E4/64 - 5*E6/256 )*Lat0-( 3*esq/8 + 3*E4/32 + 45*E6/1024)*sin(2*Lat0)+
( 15*E4/256 + 45*E6/1024 )*sin(4*Lat0)-( 35*E6/3072 )*sin(6*Lat0)
);
float M1=M0+(*Nord-FN)/K0;
float mu1=M1/(a*(1-esq/4-3*E4/64-5*E6/256));
float Lat1=mu1+(3*e1/2 - 27*pow(e1,3)/32)*sin(2*mu1)+
(21*pow(e1,2)/16 -55*pow(e1,4)/32)*sin(2*mu1)+
(151*pow(e1,3)/96)*sin(6*mu1)+(1097*pow(e1,4)/512)*sin(8*mu1);
float T1=tan(Lat1);
float C1=e_strih_sq*pow(cos(Lat1),2);
float nu1=a / sqrt(1 - esq * pow(sin(Lat1),2));
float D=(*East- FE)/(nu1*K0);
float ro1=a*(1-esq)/pow((1-esq*pow(sin(Lat1),2)),1.5);
*Nord=Lat1-(nu1*tan(Lat1)/ro1)*(pow(D,2)/2-(5+3*T1+10*C1-4*pow(C1,2)-9*e_strih_sq)*pow(D,4)/24+
(61+90*T1+298*C1+45*pow(T1,2)-252*e_strih_sq-3*pow(C1,2))*pow(D,6)/720 );
*East=Lon0+(D-(1+2*T1+C1)*pow(D,3)/6 + (5-2*C1+28*T1-3*pow(C1,2)+8*e_strih_sq+24*pow(T1,2))*pow(D,5)/120)/
cos(Lat1);
*Nord=*Nord*57.29577951-0.00037;
*East=*East*57.29577951;
}
//------------------------------------------------------------------------------
Ну и самое простое. Функция которая выдает значение Центрального меридиана зоны.
Ей даешь значение долготы, в градусах, она возвращает значение ЦМ.
int GetCM(size_t longitude)
{
int Lon0=(longitude/6)*6+3;
return Lon0;
}
Re: Не знаю может поможет..
Спасибо, но такое я уже писал давно и оно мне в данном случае не нужно...
Кто сейчас на конференции
Сейчас этот форум просматривают: нет зарегистрированных пользователей и 137 гостей