Last active
December 1, 2017 14:04
-
-
Save rmania/25150d40fd152e9adb9995d0514f74fc to your computer and use it in GitHub Desktop.
common geo manipulations
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from functools import partial | |
| import pyproj as proj | |
| from shapely.ops import transform | |
| from shapely.geometry import mapping, shape | |
| import json | |
| def rd2wgsGeojson(geojson): | |
| # convert geojson from RD new to WSG84 | |
| reprojection = partial(proj.transform, | |
| # Source coordinate system | |
| proj.Proj(init='epsg:28992'), | |
| # Destination coordinate system | |
| proj.Proj(init='epsg:4326')) | |
| g1 = shape(geojson) | |
| g2 = transform(reprojection, g1) # apply projection | |
| geom = json.dumps(mapping(g2)) | |
| return geom | |
| # -------------------------------- | |
| # same but for a Pandas series. This needs to be a string without NaN values | |
| from pyproj import Proj, transform | |
| inProj = Proj(init='epsg:28992') | |
| wgs84= Proj("+init=EPSG:4326") # LatLon with WGS84 datum used by GPS units and Google Earth | |
| for i in test[['x_coord', 'y_coord']]: | |
| test.copy().loc[:,i] = test.loc[:,i].astype(str) | |
| def convertCoords(row): | |
| x2, y2 = transform(inProj,wgs84, row['x_coord'], row['y_coord']) | |
| return pd.Series({'long':x2,'lat':y2}) | |
| # create subset with latitude/ longitude (NaN values to be removed) | |
| lon_lat = test.dropna(subset = ['x_coord', 'y_coord'])[['x_coord', 'y_coord']] | |
| lon_lat = lon_lat.join(lon_lat.apply(convertCoords, axis=1)) | |
| # join back on original frame | |
| test = pd.merge(test, lon_lat, left_index= True, right_index=True, how= 'left') | |
| ## ===================================================== | |
| # Latitude/Longitude to/from Web Mercator | |
| import math | |
| def toWGS84(xLon, yLat): | |
| # Check if coordinate out of range for Latitude/Longitude | |
| if (abs(xLon) < 180) and (abs(yLat) > 90): | |
| return | |
| # Check if coordinate out of range for Web Mercator | |
| # 20037508.3427892 is full extent of Web Mercator | |
| if (abs(xLon) > 20037508.3427892) or (abs(yLat) > 20037508.3427892): | |
| return | |
| semimajorAxis = 6378137.0 # WGS84 spheriod semimajor axis | |
| latitude = (1.5707963267948966 - (2.0 * math.atan(math.exp((-1.0 * yLat) / semimajorAxis)))) * (180/math.pi) | |
| longitude = ((xLon / semimajorAxis) * 57.295779513082323) - ((math.floor((((xLon / semimajorAxis) * 57.295779513082323) + 180.0) / 360.0)) * 360.0) | |
| return [longitude, latitude] | |
| def toWebMercator(xLon, yLat): | |
| # Check if coordinate out of range for Latitude/Longitude | |
| if (abs(xLon) > 180) and (abs(yLat) > 90): | |
| return | |
| semimajorAxis = 6378137.0 # WGS84 spheriod semimajor axis | |
| east = xLon * 0.017453292519943295 | |
| north = yLat * 0.017453292519943295 | |
| northing = 3189068.5 * math.log((1.0 + math.sin(north)) / (1.0 - math.sin(north))) | |
| easting = semimajorAxis * east | |
| return [easting, northing] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment