Create sharpmap tile like google

Topics: SharpMap v0.9 / v1.x
Developer
Sep 15, 2011 at 6:57 AM

Hi all,

I use sharp map to render tiltle image. On zoom level 0, the whole 360 degrees of longitude ( world map)  are visible in a single tile ( 256x256 )

The result the image of sharpmap and google is different.

How to render image 256x256 like google tile?. Maybe scale of map in sharpmap is different

My DPI: 96

Resolution of computer: 1024x768

View the result at: http://ada.com.vn/dmdocuments/Img_world/level0/image.html

Thanks!

TrieuVy

Coordinator
Sep 15, 2011 at 7:58 AM

You need to either

  • Reproject your data to EPSG:3857 beforehand using e.g. ogr2ogr or QGis
  • Do on the fly reprojection with SharpMap and either ProjNet or DotSpatial.Projections

Hth FObermaier

Developer
Sep 15, 2011 at 8:07 AM

Hi FObermaier,

- May be you said: EPSG: 3785 ( or EPSG:3875 )?

- I convert to EPSG: 3785 in postgis provider:

Use: INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values (
   3785, 'EPSG', 3785,
   '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ',
   'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]'
);

and

UPDATE worldcountries_boundary SET the_geom = Transform(the_geom,3785);
-- constrain
ALTER TABLE worldcountries_boundary
  ADD CONSTRAINT enforce_srid_the_geom CHECK (srid(the_geom) = 3785);
---UpdateGeometrySRID
 Select UpdateGeometrySRID('worldcountries_boundary', 'the_geom', 3785);--32648 he qui chieu thuc

- The result is the same

Please give me advice!

Coordinator
Sep 15, 2011 at 8:27 AM

yes, that is what you need to do, though the official epsg code is EPSG:3857

This has been muddeled up in the past (see: http://www.epsg-registry.org/ and query for 3785, 3857 and 3875)

The code is EPSG:3857!

cheers FObermaier

Developer
Sep 15, 2011 at 8:48 AM

:)

ok, the code is EPSG:3875

Now, I get map follow:

_map.Size = new Size(int.Parse(txtWidth.Text), int.Parse(txtHeight.Text));

_map.ZoomToExtents();
          
 Image img = _map.GetMap();

The img is the same result : http://ada.com.vn/dmdocuments/Img_world/level0/image.html

Developer
Sep 15, 2011 at 8:49 AM

_map.Size = new Size(256, 256);

_map.ZoomToExtents();
          
 Image img = _map.GetMap();

Coordinator
Sep 15, 2011 at 8:55 AM

Do you happen to know what spatial reference system the data was in before you did transform?

If it was undefined (-1) the transformation will not do anything to your data

Developer
Sep 15, 2011 at 9:24 AM

Hi

layer.CoordinateTransformation is null.

how to set CoordinateTransformation of layer is EPSG:3875?

Thanks,

TrieuVy

Coordinator
Sep 15, 2011 at 9:52 AM
Edited Sep 15, 2011 at 10:19 AM

if you are using a postgis db backend, there is no need to set the coordinate transformation, since you can do the coordinate transformation once beforehand.

Assuming you use the postgis demo db, you have for the layers countries, cities and rivers the srid 4326 set, which resembles EPSG:4326

Now if you want e.g. the countries layer in EPSG:3857 you need to do the following:

SELECT AddGeometryColumn('countries', 'geom3857', 3857, 'POLYGON', 2);
UPDATE countries SET "geom3857" = st_transform("wkb_geometry", 3857);
CREATE INDEX "idx_countries_geom3857" ON countries  USING GIST("geom3857");
SELECT DropGeometryColumn('countries', 'wkb_geometry');
--Maybe you need to drop the index as well
--DROP INDEX "idx_countries_geom";

After that you can change your layer setup to use "geom3857" geometry column

 

Hth FObermaier

Developer
Sep 15, 2011 at 10:12 AM

Thanks FObermaier!

- the first:

spatial_ref_sys is not esist 3857 srid

How to insert this srid?

TrieuVy

Coordinator
Sep 15, 2011 at 10:35 AM
Edited Sep 15, 2011 at 10:36 AM

 

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) values ( 
3857, 
'EPSG',
3857,
$L$PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]$L$,
$L$+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs$L$
);
Developer
Sep 15, 2011 at 10:50 AM

Thanks FObermaier very much!

The above insert have to change :

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, srtext,proj4text) values ( 
3857, 
'EPSG',
3857,
$L$PROJCS["Popular Visualisation CRS / Mercator (deprecated)",GEOGCS["Popular Visualisation CRS",DATUM["Popular_Visualisation_Datum",SPHEROID["Popular Visualisation Sphere",6378137,0,AUTHORITY["EPSG","7059"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6055"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4055"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],AUTHORITY["EPSG","3785"],AXIS["X",EAST],AXIS["Y",NORTH]]$L$,
$L$+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs$L$
);
May be you have a mistake between srtext and proj4text column.

After I transformed to geom3857 column
 

_map.Size = new Size(256, 256);

_map.ZoomToExtents();
          
 Image img = _map.GetMap();

The img is no change
:(
Coordinator
Sep 15, 2011 at 11:48 AM

please post the results for your data table

SELECT getsrid("<geometrycolumn>"), st_astext(st_centroid("geometrycolumn")) from "geometrytable" limit 1;

Developer
Sep 19, 2011 at 2:40 AM

Getsrid: 3857

astext: "POINT(-86.1885979305625 14.8459051091555)"

I'm using ProjNet to convert Wgs84toGoogleMercator. But some geometry transform wrong.

Coordinator
Sep 19, 2011 at 7:43 AM

then you have a mismatch in your data:

The coordinates seem to be WGS84, and the SRID says they are EPSG:3857. Therfore the approach described above won't work since postgis thinks no transformation needs to be done.

Revert your original SRID

SELECT UpdateGeometrySRID('worldcountries_boundary', 'the_geom', 4326);

After that you should perform the steps transforming to EPSG:3857:

SELECT AddGeometryColumn('countries', 'geom3857', 3857, 'POLYGON', 2);
UPDATE countries SET "geom3857" = st_transform("wkb_geometry", 3857);
CREATE INDEX "idx_countries_geom3857" ON countries  USING GIST("geom3857");
Developer
Sep 20, 2011 at 4:40 AM

I run script follow:

 

SELECT UpdateGeometrySRID('worldcountries_boundary', 'the_geom', 4326);
SELECT AddGeometryColumn('worldcountries_boundary', 'geom3857', 3857, 'MULTILINESTRING', 2);
UPDATE worldcountries_boundary SET "geom3857" = transform("the_geom", 3857);
CREATE INDEX "idx_worldcountries_boundary_geom3857" ON worldcountries_boundary  USING GIST("geom3857");

 

But some transform wrong of update sentence

UPDATE worldcountries_boundary SET "geom3857" = transform("the_geom", 3857);

-->
ERROR: transform: couldn't project point: -20 (tolerance condition error)
SQL state: XX000

Note:

worldcountries_boundary table is boundary line of worldcountries data. It is LinesString or MultiLineString
Developer
Sep 20, 2011 at 5:19 AM

There are 2 problem

- If i use worldcountries.shp in data sample of Esri .( worldcoutries in wgs84 is the same to google data)

ERROR: transform: couldn't project point: -20 (tolerance condition error)
SQL state: XX000

- If i use countries.shp in sharpmap demo data. ( countries in wgs84 is not the same to google data, difference of several meters)

Transform is ok.

The result of image 256x256 is different follow: http://ada.com.vn/dmdocuments/Img_world/level0/image.html

Coordinator
Sep 20, 2011 at 9:36 AM
Edited Sep 20, 2011 at 9:37 AM

> - If i use worldcountries.shp in data sample of Esri .( worldcoutries in wgs84 is the same to google data)

sorry, I cannot follow.

Right after you import the shapefile into postgresql/postgis, what does the geometry_tables state for the table, what is the definition for the srid in the table row?
Maybe you must not call UpdateGeometrySRID after all.

Developer
Sep 20, 2011 at 9:58 AM

This is all of scripts that i run  after i import the shapefile into postgresql/postgis

SELECT UpdateGeometrySRID('countries', 'the_geom', 4326);
update countries set the_geom = setSRID(the_geom,4326);
SELECT AddGeometryColumn('countries', 'geom3857', 3857, 'MULTIPOLYGON', 2);
UPDATE countries SET "geom3857" = transform("the_geom", 3857);
CREATE INDEX "idx_countries_geom3857" ON countries  USING GIST("geom3857");
select DropGeometryColumn('countries','the_geom'); --DropGeometryColumn(<schema_name>, <table_name>, <column_name>) 

Coordinator
Sep 20, 2011 at 10:27 AM

and what happens if you leave out the first two steps?

Developer
Sep 20, 2011 at 10:42 AM

if i leave out the first two steps. Srid curent is -1 and the script follow is error:

UPDATE countries SET "geom3857" = transform("the_geom", 3857);
--> ERROR: Input geometry has unknown (-1) SRID
SQL state: XX000

 

Coordinator
Sep 20, 2011 at 11:58 AM

Please post the contents of worldcountries_boundary.prj

Developer
Sep 21, 2011 at 2:29 AM

- worldcountries.shp  in Esri demo data, i can not find worldcountries.prj file

- countries.shp  in sharpmap demo data . Countries.prj file have the content follow:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]]

Developer
Sep 24, 2011 at 5:28 AM

Hi FObermaier!

Thanks for your support,

I used Projnet in order to convert to mecator. The result is successful!

Please view Issue:

http://sharpmap.codeplex.com/workitem/31495

Coordinator
Sep 25, 2011 at 11:48 AM

Hello Trieuvy,

I'm glad you could make it work. It looks nice :).

cheers, FObermaier