1

I have a problem/question what is the most efficient way to go from WGS84 (EPSG:4326) to Web Mercator.

Outline of the structure right now:

  1. MySQL Database: with point Latitude/Longitude Data [~200 million points] these points I bin myself with a SQL query, right now I just bin them with (I insert into separate table for convenience):

    SELECT 
        (FLOOR((longitude/100000)/(360/(POWER(2,16))))*(360/(POWER(2,16)))) AS pixel_X,     
        (FLOOR((latitude/100000)/(180/(POWER(2,15))))*(180/(POWER(2,15)))) AS pixel_Y,
        count(*) AS Amount
        FROM PointData
        GROUP BY pixel_X, pixel_Y
        ORDER BY pixel_X ASC, pixel_Y ASC;
    

    I think this should give me a raster in EPSG:4326.
    This is where I think would be the best place to introduce a change to bin into EPSG:3857 (Web Mercator, EPSG:900913). I bin the information here because it's impossible to send 200-300 million points to GDAL so I have to give up some correctness.
  2. GDAL_grid that creates a raster img-file from the data in the database (the DBConnection.vrt is just a database connection).

    gdal_grid -of HFA -a nearest:radius1=0.0054931640625:radius2=0.0054931640625 -outsize 65536 32768 -l layer DBConnection.vrt Rasters/output.img --config GDAL_NUM_THREADS ALL_CPUS


    Right now this creates an image in EPSG:4326. The nearest algorithm is used to keep the binned raster from MySQL as unchanged as possible.
  3. MapServer reading the img raster

     MAP
            NAME "TestMap"
            PROJECTION
                "init=epsg:4326"
            END
            STATUS OFF
            EXTENT -180 -90 180 90
            UNITS DD
            SHAPEPATH "./layers"
            WEB
                IMAGEPATH "/ms4w/tmp/ms_tmp/"
                IMAGEURL "/ms_tmp/"
                METADATA
                    "wms_title" "Test Server"
                    "wms_onlineresource" "http://localhost/cgi-bin/mapserv?map=/ms4w/apps/TestMaps/TestMap.map&"
                    "wms_enable_request" "*"
                END
            END         
            [...]   
            LAYER
                NAME TestLayer
                STATUS OFF
                TYPE RASTER
                DATA "Rasters\output.img"
                METADATA
                    "wms_title" "output"
                    "wms_transparent" "true"
                END
                PROJECTION
                    "init=epsg:4326"
                END
                CLASSITEM "pixel"
                PROCESSING "SCALE=0,1000"
                PROCESSING "SCALE_BUCKETS=100"
                CLASS
                    [...] (just some coloring to give a makeshift heatmap)
                END
            END
        END
        
  4. OpenLayers connects to MapServer and shows the map in the browser.

     function init() {
        map = new OpenLayers.Map("Map", {
                [...]
            ]
        });
        [...]
        var TestLayer = new OpenLayers.Layer.WMS(
            "TestLayer", 
            "http://localhost/cgi-bin/mapserv.exe?map=/ms4w/apps/TestMaps/TestMap.map", 
            {layers: 
                'TestLayer',
                transparent: 'true'
            }, {
                isBaseLayer:false,
                visibility: true
            }
        );
        [...]
    

So at this point the maps show up fine in a EPSG:4326 projection, no problem at all, I have EPSG:4326 country polygons that I use as a base layer right now.

Result in epsg:4326

The problem is I would like to have a Web Mercator map, street map and so on, but I can't get the projection to change.

So far I've tried:

  • changing the Projection "init:epsg:4326" in MapServer to epsg:3857 but to no avail, the layers are not recognized by open layers any more and I just get picture not received error.
  • using gdal_warp on the output.img and change to epsg:3857, didn't work either, same problem, also this is really time inefficient because warp takes forever and I might as well put it into the right projection right away in the database or in the GDAL_GRID.

I have probably been misunderstanding projections and GIS terms so if so let me know, this is the first time I do anything like this, only worked with the Google Maps API before with small amounts of Lat/Lng points where you just send the Lat/Lng to Google and it does the rest.

I tried to include as much information as I could, if anything is missing or unclear let me know.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.