GitHubRGI /
geopackage-python
| @@ 288-370 (lines=83) @@ | ||
| 285 | return meters_x, meters_y |
|
| 286 | ||
| 287 | ||
| 288 | class ScaledWorldMercator(EllipsoidalMercator): |
|
| 289 | """ |
|
| 290 | Scaled World Mercator projection class that holds specific calculations |
|
| 291 | and formulas for EPSG9804/9805 projection proposed by NGA Craig Rollins. |
|
| 292 | """ |
|
| 293 | ||
| 294 | def __init__(self): |
|
| 295 | """ |
|
| 296 | Constructor |
|
| 297 | """ |
|
| 298 | super(ScaledWorldMercator, self).__init__() |
|
| 299 | ||
| 300 | @staticmethod |
|
| 301 | def pixel_size(z): |
|
| 302 | """ |
|
| 303 | Calculates the pixel size for a given zoom level. |
|
| 304 | """ |
|
| 305 | return 125829.12 / 2**z |
|
| 306 | ||
| 307 | @staticmethod |
|
| 308 | def lat_to_northing(lat): |
|
| 309 | """ |
|
| 310 | Convert a latitude to a northing |
|
| 311 | / / pi phi \ / 1 - e sin(phi) \ e/2 \ |
|
| 312 | y(phi) = R ln| tan| --- + --- | | -------------- | | |
|
| 313 | \ \ 4 2 / \ 1 + e sin(phi) / / |
|
| 314 | """ |
|
| 315 | r = 6378137.0 * 0.857385503731176 |
|
| 316 | e = 0.081819190842621 |
|
| 317 | return r * log(tan((pi / 2 + lat) / 2) * ((1 - e * sin(lat)) / |
|
| 318 | (1 + e * sin(lat)))**(e / 2)) |
|
| 319 | ||
| 320 | @staticmethod |
|
| 321 | def tile_to_lat_lon(z, x, y): |
|
| 322 | """ |
|
| 323 | Returns the lat/lon coordinates of the bottom-left corner of the input |
|
| 324 | tile. Finds the value numerically (using the secant method). A scale |
|
| 325 | factor has been added specifically for scaled world mercator. |
|
| 326 | ||
| 327 | Inputs: |
|
| 328 | z -- zoom level value for input tile |
|
| 329 | x -- tile column value for input tile |
|
| 330 | y -- tile row value for input tile |
|
| 331 | """ |
|
| 332 | n = 2.0**z |
|
| 333 | r = 6378137.0 * 0.857385503731176 |
|
| 334 | lon = x / n * 360.0 - 180.0 |
|
| 335 | my = (y - 2**(z - 1)) * r * pi * 2 / 2**z |
|
| 336 | ||
| 337 | def f(phi): |
|
| 338 | return ScaledWorldMercator.lat_to_northing(phi) - my |
|
| 339 | ||
| 340 | lat = 0.0 |
|
| 341 | oldLat = 1.0 |
|
| 342 | diff = 1.0 |
|
| 343 | while abs(diff) > 0.0001: |
|
| 344 | newLat = lat - f(lat) * (lat - oldLat) / (f(lat) - f(oldLat)) |
|
| 345 | if newLat > 1.4849969713855238: |
|
| 346 | newLat = 1.4849969713855238 |
|
| 347 | elif newLat < -1.4849969713855238: |
|
| 348 | newLat = -1.4849969713855238 |
|
| 349 | oldLat = lat |
|
| 350 | lat = newLat |
|
| 351 | diff = lat - oldLat |
|
| 352 | lat = lat * 180.0 / pi |
|
| 353 | return lat, lon |
|
| 354 | ||
| 355 | def tile_to_meters(self, z, x, y): |
|
| 356 | """ |
|
| 357 | Returns the meter coordinates of the bottom-left corner of the input |
|
| 358 | tile. A scale factor has been added to the longitude meters |
|
| 359 | calculation. |
|
| 360 | ||
| 361 | Inputs: |
|
| 362 | z -- zoom level value for input tile |
|
| 363 | x -- tile column (longitude) value for input tile |
|
| 364 | y -- tile row (latitude) value for input tile |
|
| 365 | """ |
|
| 366 | lat, lon = self.tile_to_lat_lon(z, x, y) |
|
| 367 | meters_x = lon * (pi * (6378137.0 * 0.857385503731176)) / 180.0 |
|
| 368 | meters_y = self.lat_to_northing(lat * pi / 180.0) |
|
| 369 | # Instituting a 2 decimal place round to ensure accuracy |
|
| 370 | return meters_x, round(meters_y, 2) |
|
| 371 | ||
| 372 | ||
| 373 | class ZoomMetadata(object): |
|
| @@ 214-285 (lines=72) @@ | ||
| 211 | return '%.7f' % (int(coord * 10000000) / float(10000000)) |
|
| 212 | ||
| 213 | ||
| 214 | class EllipsoidalMercator(Mercator): |
|
| 215 | """ |
|
| 216 | Ellipsoidal Mercator projection class that holds specific calculations and |
|
| 217 | formulas for EPSG3395. |
|
| 218 | """ |
|
| 219 | ||
| 220 | def __init__(self): |
|
| 221 | """ |
|
| 222 | Constructor |
|
| 223 | """ |
|
| 224 | super(EllipsoidalMercator, self).__init__() |
|
| 225 | ||
| 226 | @staticmethod |
|
| 227 | def lat_to_northing(lat): |
|
| 228 | """ |
|
| 229 | Convert a latitude to a northing |
|
| 230 | / / pi phi \ / 1 - e sin(phi) \ e/2 \ |
|
| 231 | y(phi) = R ln| tan| --- + --- | | -------------- | | |
|
| 232 | \ \ 4 2 / \ 1 + e sin(phi) / / |
|
| 233 | """ |
|
| 234 | r = 6378137.0 |
|
| 235 | e = 0.081819190842621 |
|
| 236 | return r * log(tan((pi / 2 + lat) / 2) * ((1 - e * sin(lat)) / |
|
| 237 | (1 + e * sin(lat)))**(e / 2)) |
|
| 238 | ||
| 239 | @staticmethod |
|
| 240 | def tile_to_lat_lon(z, x, y): |
|
| 241 | """ |
|
| 242 | Returns the lat/lon coordinates of the bottom-left corner of the input |
|
| 243 | tile. Finds the value numerically (using the secant method). |
|
| 244 | ||
| 245 | Inputs: |
|
| 246 | z -- zoom level value for input tile |
|
| 247 | x -- tile column value for input tile |
|
| 248 | y -- tile row value for input tile |
|
| 249 | """ |
|
| 250 | n = 2.0**z |
|
| 251 | lon = x / n * 360.0 - 180.0 |
|
| 252 | my = (y - 2**(z - 1)) * 6378137 * pi * 2 / 2**z |
|
| 253 | ||
| 254 | def f(phi): |
|
| 255 | return EllipsoidalMercator.lat_to_northing(phi) - my |
|
| 256 | ||
| 257 | lat = 0.0 |
|
| 258 | oldLat = 1.0 |
|
| 259 | diff = 1.0 |
|
| 260 | while abs(diff) > 0.0001: |
|
| 261 | newLat = lat - f(lat) * (lat - oldLat) / (f(lat) - f(oldLat)) |
|
| 262 | if newLat > 1.48499697138: |
|
| 263 | newLat = 1.48499697138 |
|
| 264 | elif newLat < -1.48499697138: |
|
| 265 | newLat = -1.48499697138 |
|
| 266 | oldLat = lat |
|
| 267 | lat = newLat |
|
| 268 | diff = lat - oldLat |
|
| 269 | lat = lat * 180.0 / pi |
|
| 270 | return lat, lon |
|
| 271 | ||
| 272 | def tile_to_meters(self, z, x, y): |
|
| 273 | """ |
|
| 274 | Returns the meter coordinates of the bottom-left corner of the input |
|
| 275 | tile. |
|
| 276 | ||
| 277 | Inputs: |
|
| 278 | z -- zoom level value for input tile |
|
| 279 | x -- tile column (longitude) value for input tile |
|
| 280 | y -- tile row (latitude) value for input tile |
|
| 281 | """ |
|
| 282 | lat, lon = self.tile_to_lat_lon(z, x, y) |
|
| 283 | meters_x = lon * self.origin_shift / 180.0 |
|
| 284 | meters_y = self.lat_to_northing(lat * pi / 180.0) |
|
| 285 | return meters_x, meters_y |
|
| 286 | ||
| 287 | ||
| 288 | class ScaledWorldMercator(EllipsoidalMercator): |
|
