||
问题:我国公开使用的在线地图都是加偏(火星坐标GCJ02)的,所以想要正确显示矢量shapefile文件,就必须进行加偏显示。
解决:在MapWinGIS中增加WGS84转GCJ02功能;
一、在项目中新建GCJ02筛选器,然后把gcj02.h和gcj02.cpp添加进去
//gcj02.h文件
#ifndef TRANSFORM_HEADER
#define TRANSFORM_HEADER
void wgs2gcj(double wgsLat, double wgsLng, double *gcjLat, double *gcjLng);
void gcj2wgs(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng);
void gcj2wgs_exact(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng);
double distance(double latA, double lngA, double latB, double lngB);
#endif
//gcj02.cpp文件
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#include "gcj02.h"
#include "stdafx.h"
#undef INLINE
#if defined(__STDC_VERSION__) && __STDC_VERSION__ >= 199900L
#define INLINE inline
#else
#define INLINE
#endif /* STDC */
#define fabs(x) __ev_fabs(x)
/* do not use >=, or compiler may not know to optimize this with a simple (x & 0x80...0). */
INLINE static double __ev_fabs(double x){ return x > 0.0 ? x : -x; }
INLINE static int outOfChina(double lat, double lng) {
if (lng < 72.004 || lng > 137.8347 || lat < 0.8293 || lat > 55.8271) {
return 1;
}
return 0;
}
#define EARTH_R 6378137.0
void xy_latlng_transform(double x, double y, double *lat, double *lng) {
double xy = x * y;
double absX = sqrt(fabs(x));
double xPi = x * M_PI;
double yPi = y * M_PI;
double d = 20.0*sin(6.0*xPi) + 20.0*sin(2.0*xPi);
*lat = d;
*lng = d;
*lat += 20.0*sin(yPi) + 40.0*sin(yPi/3.0);
*lng += 20.0*sin(xPi) + 40.0*sin(xPi/3.0);
*lat += 160.0*sin(yPi/12.0) + 320*sin(yPi/30.0);
*lng += 150.0*sin(xPi/12.0) + 300.0*sin(xPi/30.0);
*lat *= 2.0 / 3.0;
*lng *= 2.0 / 3.0;
*lat += -100.0 + 2.0*x + 3.0*y + 0.2*y*y + 0.1*xy + 0.2*absX;
*lng += 300.0 + x + 2.0*y + 0.1*x*x + 0.1*xy + 0.1*absX;
}
static void delta(double lat, double lng, double *dLat, double *dLng) {
if ((dLat == NULL) || (dLng == NULL)) {
return;
}
const double ee = 0.00669342162296594323;
xy_latlng_transform(lng - 105.0, lat - 35.0, dLat, dLng);
double radLat = lat / 180.0 * M_PI;
double magic = sin(radLat);
magic = 1 - ee*magic*magic;
double sqrtMagic = sqrt(magic);
*dLat = (*dLat * 180.0) / ((EARTH_R * (1 - ee)) / (magic * sqrtMagic) * M_PI);
*dLng = (*dLng * 180.0) / (EARTH_R / sqrtMagic * cos(radLat) * M_PI);
}
void wgs2gcj(double wgsLat, double wgsLng, double *gcjLat, double *gcjLng) {
if ((gcjLat == NULL) || (gcjLng == NULL)) {
return;
}
if (outOfChina(wgsLat, wgsLng)) {
*gcjLat = wgsLat;
*gcjLng = wgsLng;
return;
}
double dLat, dLng;
delta(wgsLat, wgsLng, &dLat, &dLng);
*gcjLat = wgsLat + dLat;
*gcjLng = wgsLng + dLng;
}
void gcj2wgs(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng) {
if ((wgsLat == NULL) || (wgsLng == NULL)) {
return;
}
if (outOfChina(gcjLat, gcjLng)) {
*wgsLat = gcjLat;
*wgsLng = gcjLng;
return;
}
double dLat, dLng;
delta(gcjLat, gcjLng, &dLat, &dLng);
*wgsLat = gcjLat - dLat;
*wgsLng = gcjLng - dLng;
}
void gcj2wgs_exact(double gcjLat, double gcjLng, double *wgsLat, double *wgsLng) {
double initDelta = 0.01;
double threshold = 0.000001;
double dLat, dLng, mLat, mLng, pLat, pLng;
dLat = dLng = initDelta;
mLat = gcjLat - dLat;
mLng = gcjLng - dLng;
pLat = gcjLat + dLat;
pLng = gcjLng + dLng;
int i;
for (i=0; i<30; i++) {
*wgsLat = (mLat+pLat) / 2;
*wgsLng = (mLng+pLng) / 2;
double tmpLat, tmpLng;
wgs2gcj(*wgsLat, *wgsLng, &tmpLat, &tmpLng);
dLat = tmpLat - gcjLat;
dLng = tmpLng - gcjLng;
if ( (fabs(dLat)<threshold) && (fabs(dLng)<threshold) ) {
return;
}
if (dLat > 0) {
pLat = *wgsLat;
} else {
mLat = *wgsLat;
}
if (dLng > 0) {
pLng = *wgsLng;
} else {
mLng = *wgsLng;
}
}
}
二、在Shapefile.cpp文件中修改如下:
1)增加头文件
#include "gcj02.h"
2)ReprojectCore()函数中增加过渡坐标系统WGS84
//--------------------chenchao-------------------
OGRSpatialReference projTemp;
projTemp.importFromEPSG(4326);//WGS84
OGRCoordinateTransformation* transfST = OGRCreateCoordinateTransformation(projSource, &projTemp);
OGRCoordinateTransformation* transf = OGRCreateCoordinateTransformation(&projTemp, projTarget);
//--------------------chenchao-------------------
//--------------------chenchao-------------------
BOOL rest = transfST->Transform(numPoints, x, y);
double lon, lat;
for (long j = 0; j < numPoints; j++)
{
wgs2gcj(y[j], x[j], &lat, &lon);
x[j] = lon;
y[j] = lat;
}
BOOL res = transf->Transform( numPoints, x, y);
//--------------------chenchao-------------------
三、转换前后的效果
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-23 19:05
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社