Calculate row and column for WMTS request from Lat/Long

I'm writing a request to a WMTS server that has an upper left corner at 90, -180 and spans the whole Earth. I want to get a tile containing some point (targetLat, targetLong).

Suppose we have some tile matrix with a given matrixWidth, matrixHeight, tileWidth, and tileHeight, and some upper-left corner matrxLat, matrixLon.

Naively, one could calculate the N-S distance from the upper-left corner of the matrix to the target point by taking the haversine distance from (targetLat, targetLong) to (matrixLat, targetLong), and the E-W distance from (targetLat, targetLong) to (targetLat, matrixLong), then use simple algebra to find the correct tile from there.

The problem is that the E-W distance is going to change depending on which latitude I measure it from, so I don't think this will work. How can I translate my target lat-long point into a meter offset from the matrix corner in the correct coordinates?

• It seems that you have been reading the standard portal.opengeospatial.org/files/?artifact_id=35326 and look at the image about the tile space on page 24. That's fine. What you need to do is to convert your lat, long point to a point that is in the same coordinate system than the tile matrix set that you are going to use. If I read right you want to use a matrix set with metric units. When you have your reference point in meters as well you can start calculating by the upper left corner, pixels per tile, and meters per pixel. – user30184 Oct 30 '18 at 7:16
• @user30184 Converting my lat, long to the correct coordinate system is, I think my entire problem. I tried to edit to make things more clear. – JETM Oct 30 '18 at 11:20
• It might be easiest to use some existing coordinate transformation library like Proj4 proj4.org, Proj4js proj4js.org or GeoTools. – user30184 Oct 30 '18 at 11:48

My issue came from assuming that each pixel mapped to a uniform physical distance, and therefore each tile has the same dimensions measured in meters. This is of course absurd.

Instead, if we have EPSG:4326 coordinate system, we can assume that each tile has the same dimensions as measured in degrees latitude and longitude.

Given an upper left corner (-180, 90), matrix width of MX tiles, matrix height of MY tiles, and point (Y, X), we have

tiles_per_degree_longitude = MX / 360
tiles_per_degree_latitude = MY / 180
dX = X + 180
dY = 90 - Y
tileX = dX * tiles_per_degree_longitude
tileY = dY * tiles_per_degree_latitude

I have tested the tiles returned from this calculation against the corresponding latitude/longitude on Google Maps.

I assume this does not work if you do not have an EPSG:4326 coordinate system, but I do not know for sure.