Commit 6c2f6c1e authored by Jon Pye's avatar Jon Pye
Browse files

Committing initial play files

parents
"""
The effect of badly referencing an ellipse
------------------------------------------
This example demonstrates the effect of referencing your data to an incorrect
ellipse.
First we define two coordinate systems - one using the World Geodetic System
established in 1984 and the other using a spherical globe. Next we extract
data from the Natural Earth land dataset and convert the Geodetic
coordinates (referenced in WGS84) into the respective coordinate systems
that we have defined. Finally, we plot these datasets onto a map assuming
that they are both referenced to the WGS84 ellipse and compare how the
coastlines are shifted as a result of referencing the incorrect ellipse.
"""
import cartopy.crs as ccrs
import cartopy.feature
from cartopy.io.img_tiles import MapQuestOpenAerial
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D as Line
from matplotlib.patheffects import Stroke
import numpy as np
import shapely.geometry as sgeom
from shapely.ops import transform as geom_transform
def transform_fn_factory(target_crs, source_crs):
"""
Return a function which can be used by ``shapely.op.transform``
to transform the coordinate points of a geometry.
The function explicitly *does not* do any interpolation or clever
transformation of the coordinate points, so there is no guarantee
that the resulting geometry would make any sense.
"""
def transform_fn(x, y, z=None):
new_coords = target_crs.transform_points(source_crs,
np.asanyarray(x),
np.asanyarray(y))
return new_coords[:, 0], new_coords[:, 1], new_coords[:, 2]
return transform_fn
def main():
# Define the two coordinate systems with different ellipses.
wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
ellipse='WGS84'))
sphere = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
ellipse='sphere'))
# Define the coordinate system of the data we have from Natural Earth and
# acquire the 1:10m physical coastline shapefile.
geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
dataset = cartopy.feature.NaturalEarthFeature(category='physical',
name='coastline',
scale='10m')
# Create a MapQuest map tiler instance, and use its CRS for the GeoAxes.
tiler = MapQuestOpenAerial()
ax = plt.axes(projection=tiler.crs)
plt.title('The effect of incorrectly referencing the Solomon Islands')
# Pick the area of interest. In our case, roughly the Solomon Islands, and
# get hold of the coastlines for that area.
extent = (155, 163, -11.5, -6)
ax.set_extent(extent, geodetic)
geoms = list(dataset.intersecting_geometries(extent))
# Add the MapQuest aerial imagery at zoom level 7.
ax.add_image(tiler, 7)
# Transform the geodetic coordinates of the coastlines into the two
# projections of differing ellipses.
wgs84_geoms = [geom_transform(transform_fn_factory(wgs84, geodetic),
geom) for geom in geoms]
sphere_geoms = [geom_transform(transform_fn_factory(sphere, geodetic),
geom) for geom in geoms]
# Using these differently referenced geometries, assume that they are
# both referenced to WGS84.
ax.add_geometries(wgs84_geoms, wgs84, edgecolor='white', color='none')
ax.add_geometries(sphere_geoms, wgs84, edgecolor='gray', color='none')
# Create a legend for the coastlines.
legend_artists = [Line([0], [0], color=color, linewidth=3)
for color in ('white', 'gray')]
legend_texts = ['Correct ellipse\n(WGS84)', 'Incorrect ellipse\n(sphere)']
legend = plt.legend(legend_artists, legend_texts, fancybox=True,
loc='lower left', framealpha=0.75)
legend.legendPatch.set_facecolor('wheat')
# Create an inset GeoAxes showing the location of the Solomon Islands.
sub_ax = plt.axes([0.7, 0.625, 0.2, 0.2], projection=ccrs.PlateCarree())
sub_ax.set_extent([110, 180, -50, 10], geodetic)
# Make a nice border around the inset axes.
effect = Stroke(linewidth=4, foreground='wheat', alpha=0.5)
sub_ax.outline_patch.set_path_effects([effect])
# Add the land, coastlines and the extent of the Solomon Islands.
sub_ax.add_feature(cartopy.feature.LAND)
sub_ax.coastlines()
extent_box = sgeom.box(extent[0], extent[2], extent[1], extent[3])
sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), color='none',
edgecolor='blue', linewidth=2)
plt.show()
if __name__ == '__main__':
main()
\ No newline at end of file
# Imports to get the data out of GeoServer
from owslib.wfs import WebFeatureService
import geojson
import numpy as np
# Imports to plot the map
import cartopy.crs as ccrs
from cartopy.io.img_tiles import MapQuestOpenAerial, MapQuestOSM, OSM
import matplotlib.pyplot as plt
from shapely.geometry import Polygon, mapping, asShape, shape, asPoint, point
def flip_geojson_coordinates(geo):
if isinstance(geo, dict):
for k, v in geo.iteritems():
if k == "coordinates":
z = np.asarray(geo[k])
f = z.flatten()
geo[k] = np.dstack((f[1::2], f[::2])).reshape(z.shape).tolist()
else:
flip_geojson_coordinates(v)
elif isinstance(geo, list):
for k in geo:
flip_geojson_coordinates(k)
wfs = WebFeatureService('http://global.devl.oceantrack.org:8080/geoserver/wfs?', version="2.0.0")
print list(wfs.contents)
a = wfs.contents['otn:animals']
b = a.boundingBoxWGS84
json_response = wfs.getfeature(typename=['otn:animals'], propertyname=None, outputFormat='application/json').read()
geom = geojson.loads(json_response)
flip_geojson_coordinates(geom)
print geom.keys()
print geom['type'], geom['crs'] # check our coordinate reference system
geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
# Open Source Imagery from MapQuest (max zoom = 16?)
tiler = MapQuestOpenAerial()
# Open Street Map (max zoom = 18?)
tiler = OSM()
ax = plt.axes(projection=tiler.crs)
dx=b[2]-b[0]
dy=b[3]-b[1]
extent = (b[0]-0.1*dx,b[2]+0.1*dx,b[1]-0.1*dy,b[3]+0.1*dy)
ax.set_extent(extent, geodetic)
ax.add_image(tiler, 7)
#ax.add_geometries([polygon],ccrs.PlateCarree(),
# facecolor=BLUE, edgecolor=GRAY,alpha=0.5)
for p in geom.get("features", []):
# How to feed this to the ax.?
#multi_poly = asShape(p.get("geometry"))
thePoint = asPoint(p.get("geometry"))
ax.add_geometries(p.get("geometry"),ccrs.PlateCarree(),
edgecolor='black',facecolor='none', color='black')
#title(name)
print 'done with geometry'
gl=ax.gridlines(draw_labels=True)
gl.xlabels_top = False
gl.ylabels_right = False
#ax.add_feature(coast_10m,edgecolor='black')
#ax.coastlines()
plt.show()
Cartopy==0.12.x
Cython==0.20.2
OWSLib==0.8.8
Pillow==2.5.1
PySAL==1.7.0
Shapely==1.3.2
backports.ssl-match-hostname==3.4.0.2
certifi==14.05.14
descartes==1.0.1
geojson==1.0.7
gsconfig==0.6.9
httplib2==0.9
matplotlib==1.3.1
nose==1.3.3
numpy==1.8.1
pandas==0.14.1
psycopg2==2.5.3
pyparsing==2.0.2
pyshp==1.2.1
python-dateutil==2.2
pytz==2014.4
scipy==0.14.0
six==1.7.3
tornado==4.0
wsgiref==0.1.2
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment