chaoshu的个人博客分享 http://blog.sciencenet.cn/u/chaoshu

博文

MapWinGIS实现GCJ02火星坐标转换

已有 3344 次阅读 2019-8-31 16:34 |个人分类:GIS|系统分类:科研笔记

问题:我国公开使用的在线地图都是加偏(火星坐标GCJ02)的,所以想要正确显示矢量shapefile文件,就必须进行加偏显示。

解决:在MapWinGIS中增加WGS84转GCJ02功能;

一、在项目中新建GCJ02筛选器,然后把gcj02.h和gcj02.cpp添加进去

image.png

image.png

//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

image.png


//--------------------chenchao-------------------

OGRSpatialReference projTemp;

projTemp.importFromEPSG(4326);//WGS84

OGRCoordinateTransformation* transfST = OGRCreateCoordinateTransformation(projSource, &projTemp);

OGRCoordinateTransformation* transf = OGRCreateCoordinateTransformation(&projTemp, projTarget);

//--------------------chenchao-------------------


image.png

//--------------------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-------------------


三、转换前后的效果

image.png

image.png



https://blog.sciencenet.cn/blog-858128-1196050.html

上一篇:MapWinGIS加载地图缓慢问题解决
下一篇:基于MFC的MapWinGIS加载tif/shapefile文件
收藏 IP: 220.180.184.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-16 19:09

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部