Тема

[X-Cray]

Киньте ссылкой на теорию, плизз :)
Есть N точек с растровыми координатами. Это N ооочень велико :)
Знаю датум, знаю проекцию, знаю геодезические координаты нескольки таких точек.
Нужно найти такие координаты для всех остальных. Тоесть сделать то, что умеет Ози после привязки растра. Насколько я понял, нужно по известным данным посчитать параметры аффинного преобразования, а потом его выполнить. И что-то еще...

[BrainDrain]

Читай тут http://www-ec.njit.edu/surveying/Docume ... .pdf но если задача однократная, то проще сделать это какой-нибудь ГИС, это многие умеют (рекомендую GeoMedia, ArcGIS)

[X-Cray]

Задача однократная, но растра как такового нету, есть просто куча координат в базе, которые можно считать растровыми. Конечно, можно по ним сформировать картинку. Но, думаю, удобнее будет вывести формулы преобразования и прогнать через них все мои точки, тем более, что их сразу же надо положить в другую базу.

[X-Cray]

Спасибо за ссылку. Полезная информация

[Iitta]

Если умеете немного программировать, то задача довольно простая.
Пишете линейные уравнения (или квадратичные, если есть искажения), связывающие координаты растра и топо (типа: 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-Cray]

С программированием не расстаюсь уже лет 6 :). Вспомнилась лаба по интерполяции, что мы делали года 3 назад. Кстати, для преобразования растровых координат в координаты проекции Меркатора на сферу нужна нелинейная интерполяция?

[Iitta]

Если точки покрывают континенты, то лучше ждать ответа от монстров форума. На маленьких площадях вся меркаторность и эллипсоидальность ИМХО уже учтена в координатах известных точек. Впрочем, даже для континентов, если точки с известными координатами более или менее равномерно покрывают интересующую площадь, можно разбить всю область на небольшие куски и выполнить линейное преобразование внутри каждого из них. Т.е. прокластеризовать пространство вокруг известных точек.

[Olexa]

: Есть N точек с растровыми координатами. Это N ооочень велико :)
: Знаю датум, знаю проекцию, знаю геодезические координаты нескольки таких точек.
: Нужно найти такие координаты для всех остальных. Тоесть сделать то, что умеет Ози после привязки растра. Насколько я понял, нужно по известным данным посчитать параметры аффинного преобразования, а потом его выполнить. И что-то еще...

1. Пересчитать географические координаты точек привязки в датум и координаты проекции (которые известны).
2. Рассчитать параметры афинного преобразования по известным парам координат ("растровых" и координат проекции) для точек привязки.
3. При помощи афинного преобразования рассчитать координаты проекции для остальных точек.
4. Пересчитать координаты проеции для остальных точек в географические координаты и соответствующий датум.

[gorbva]

Вот надысь сделал аддон к мапедиту где преобразовывал координаты из 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;
}

[X-Cray]

Спасибо, но такое я уже писал давно и оно мне в данном случае не нужно...