source: subversion/applications/utils/export/osm2pgsql-intarray/reprojection.c @ 28719

Last change on this file since 28719 was 25075, checked in by jonb, 9 years ago

osm2pgsql: Clip latitude of nodes to +-85.07 degrees. This prevents the South pole from showing up at (0,0). Patch originally from Komzpa. closes #3452

File size: 6.0 KB
Line 
1/* reprojection.c
2 *
3 * Convert OSM coordinates to another coordinate system for
4 * the database (usually convert lat/lon to Spherical Mercator
5 * so Mapnik doesn't have to).
6 */
7
8#include <stdio.h>
9#include <string.h>
10#include <unistd.h>
11#include <proj_api.h>
12#include <math.h>
13
14#include "reprojection.h"
15
16/** must match expire.tiles.c */
17#define EARTH_CIRCUMFERENCE              40075016.68
18
19/** The projection of the source data. Always lat/lon (EPSG:4326). */
20static projPJ pj_source = NULL;
21
22/** The target projection (used in the PostGIS tables). Controlled by the -l/-M/-m/-E options. */
23static projPJ pj_target = NULL;
24
25/** The projection used for tiles. Currently this is fixed to be Spherical
26 *  Mercator. You will usually have tiles in the same projection as used
27 *  for PostGIS, but it is theoretically possible to have your PostGIS data
28 *  in, say, lat/lon but still create tiles in Spherical Mercator.
29 */
30static projPJ pj_tile = NULL;
31
32static int Proj;
33
34const struct Projection_Info Projection_Info[] = {
35  [PROJ_LATLONG] = { 
36     .descr     = "Latlong", 
37     .proj4text = "+init=epsg:4326",
38     .srs       = 4326, 
39     .option    = "-l" },
40  [PROJ_MERC] = { 
41     .descr     = "WGS84 Mercator", 
42     .proj4text = "+proj=merc +datum=WGS84  +k=1.0 +units=m +over +no_defs", 
43     .srs       = 3395, 
44     .option    = "-M" },
45  [PROJ_SPHERE_MERC] = { 
46     .descr     = "Spherical Mercator", 
47     .proj4text = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs", 
48     .srs       = 900913, 
49     .option    = "-m" }
50};
51static struct Projection_Info custom_projection;
52
53/** defined in expire-tiles.c; depends on the zoom level selected for expiry. */
54extern int map_width; 
55
56// Positive numbers refer the to the table above, negative numbers are
57// assumed to refer to EPSG codes and it uses the proj4 to find those.
58void project_init(int proj)
59{
60    char buffer[32];
61    Proj = proj;
62   
63    /* hard-code the source projection to be lat/lon, since OSM XML always
64     * has coordinates in degrees. */
65    pj_source = pj_init_plus("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs");
66
67    /* hard-code the tile projection to be spherical mercator always.
68     * theoretically this could be made selectable but not all projections
69     * lend themselves well to making tiles; non-spherical mercator tiles
70     * are uncharted waters in OSM. */
71    pj_tile = pj_init_plus(Projection_Info[PROJ_SPHERE_MERC].proj4text);
72
73    /* now set the target projection - the only one which is really variable */
74    if (proj >= 0 && proj < PROJ_COUNT)
75    {
76        pj_target = pj_init_plus(Projection_Info[proj].proj4text);
77    }
78    else if (proj < 0)
79    {
80        if (snprintf(buffer, sizeof(buffer), "+init=epsg:%d", -proj ) >= (int)sizeof(buffer))
81        {
82            fprintf(stderr, "Buffer overflow computing proj4 initialisation string\n");
83            exit(1);
84        }
85        pj_target = pj_init_plus(buffer);
86        if (!pj_target)
87        {
88            fprintf (stderr, "Couldn't read EPSG definition (do you have /usr/share/proj/epsg?)\n");
89            exit(1);
90        }
91    }
92           
93    if (!pj_source || !pj_target || !pj_tile) 
94    {
95        fprintf(stderr, "Projection code failed to initialise\n");
96        exit(1);
97    }
98   
99    if (proj >= 0)
100        return;
101
102    custom_projection.srs = -proj;
103    custom_projection.proj4text = pj_get_def(pj_target, 0);
104    if (snprintf(buffer, sizeof(buffer), "EPSG:%d", -proj) >= (int)sizeof(buffer))
105    {
106        fprintf(stderr, "Buffer overflow computing projection description\n");
107        exit(1);
108    }
109    custom_projection.descr = strdup(buffer);
110    custom_projection.option = "-E";
111    return;
112}
113
114void project_exit(void)
115{
116    pj_free(pj_source);
117    pj_source = NULL;
118    pj_free(pj_target);
119    pj_target = NULL;
120}
121
122struct Projection_Info const *project_getprojinfo(void)
123{
124  if( Proj >= 0 )
125    return &Projection_Info[Proj];
126  else
127    return &custom_projection;
128}
129
130void reproject(double *lat, double *lon)
131{
132    double x[1], y[1], z[1];
133   
134    /** Caution: This section is only correct if the source projection is lat/lon;
135     *  so even if it looks like pj_source was just a variable, things break if
136     *  pj_source is something else than lat/lon. */
137
138    if (Proj == PROJ_LATLONG)
139        return;
140
141    if (Proj == PROJ_SPHERE_MERC)
142    {
143      /* The latitude co-ordinate is clipped at slightly larger than the 900913 'world'
144       * extent of +-85.0511 degrees to ensure that the points appear just outside the
145       * edge of the map. */
146
147        if (*lat > 85.07)
148            *lat = 85.07;
149        if (*lat < -85.07)
150            *lat = -85.07;
151
152        *lat = log(tan(M_PI/4.0 + (*lat) * DEG_TO_RAD / 2.0)) * EARTH_CIRCUMFERENCE/(M_PI*2);
153        *lon = (*lon) * EARTH_CIRCUMFERENCE / 360.0;
154        return;
155    }
156
157    x[0] = *lon * DEG_TO_RAD;
158    y[0] = *lat * DEG_TO_RAD;
159    z[0] = 0;
160
161    /** end of "caution" section. */
162   
163    pj_transform(pj_source, pj_target, 1, 1, x, y, z);
164   
165    //printf("%.4f\t%.4f -> %.4f\t%.4f\n", *lat, *lon, y[0], x[0]);
166    *lat = y[0];
167    *lon = x[0];
168}
169
170/**
171 * Converts from (target) coordinates to tile coordinates.
172 *
173 * The zoom level for the coordinates is implicitly given in the global
174 * variable map_width.
175 */
176void coords_to_tile(double *tilex, double *tiley, double lon, double lat)
177{
178    double x[1], y[1], z[1];
179   
180    x[0] = lon;
181    y[0] = lat;
182    z[0] = 0;
183
184    if (Proj == PROJ_LATLONG)
185    {
186        x[0] *= DEG_TO_RAD;
187        y[0] *= DEG_TO_RAD;
188    }
189
190    /* since pj_tile is always spherical merc, don't bother doing anything if
191     *  destination proj is the same. */
192
193    if (Proj != PROJ_SPHERE_MERC)
194    {
195        pj_transform(pj_target, pj_tile, 1, 1, x, y, z);
196        /** FIXME: pj_transform could fail if coordinates are outside +/- 85 degrees latitude */
197    }
198   
199    /* if ever pj_tile were allowed to be PROJ_LATLONG then results would have to
200     *  be divided by DEG_TO_RAD here. */
201
202    *tilex = map_width * (0.5 + x[0] / EARTH_CIRCUMFERENCE);
203    *tiley = map_width * (0.5 - y[0] / EARTH_CIRCUMFERENCE);
204}
205
Note: See TracBrowser for help on using the repository browser.