mirror of
https://github.com/wassname/PyCRS.git
synced 2026-06-27 16:10:29 +08:00
Karim Bahgat
380 lines
11 KiB
Python
380 lines
11 KiB
Python
import pycrs
|
|
import traceback
|
|
import logging
|
|
|
|
|
|
|
|
###########################
|
|
# Drawing routine for testing
|
|
raw = None
|
|
def render_world(crs, savename):
|
|
import urllib2
|
|
import json
|
|
import pygeoj
|
|
import pyagg
|
|
import pyproj
|
|
import random
|
|
|
|
# load world borders
|
|
global raw
|
|
if not raw:
|
|
raw = urllib2.urlopen("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json").read()
|
|
rawdict = json.loads(raw)
|
|
data = pygeoj.load(data=rawdict)
|
|
|
|
# convert coordinates
|
|
fromproj = pyproj.Proj("+init=EPSG:4326")
|
|
toproj = pyproj.Proj(crs.to_proj4())
|
|
for feat in data:
|
|
if feat.geometry.type == "Polygon":
|
|
feat.geometry.coordinates = [zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1]))
|
|
for ring in feat.geometry.coordinates]
|
|
elif feat.geometry.type == "MultiPolygon":
|
|
feat.geometry.coordinates = [
|
|
[zip(*pyproj.transform(fromproj, toproj, zip(*ring)[0], zip(*ring)[1]))
|
|
for ring in poly]
|
|
for poly in feat.geometry.coordinates]
|
|
feat.geometry.update_bbox() # important to clear away old bbox
|
|
# get zoom area
|
|
data.add_all_bboxes()
|
|
data.update_bbox()
|
|
bbox = data.bbox
|
|
|
|
## # to avoid inf bounds and no render in satellite view
|
|
## xmins, ymins, xmaxs, ymaxs = zip(*(feat.geometry.bbox for feat in data))
|
|
## inf = float("inf")
|
|
## xmaxs = (xmax for xmax in xmaxs if xmax != inf)
|
|
## ymaxs = (ymax for ymax in ymaxs if ymax != inf)
|
|
## bbox = (min(xmins), min(ymins), max(xmaxs), max(ymaxs))
|
|
|
|
# set up drawing
|
|
c = pyagg.Canvas(1000,1000)
|
|
c.geographic_space()
|
|
c.zoom_bbox(*bbox)
|
|
c.zoom_out(1.3)
|
|
|
|
# draw countries
|
|
for feat in data:
|
|
try: c.draw_geojson(feat.geometry,
|
|
fillcolor=tuple(random.randrange(255) for _ in range(3)),
|
|
outlinecolor="white")
|
|
except:
|
|
# NOTE: feat.__geo_interface__ is one level too high maybe??
|
|
print("unable to draw?", feat.geometry)
|
|
|
|
# draw text of the proj4 string used
|
|
c.percent_space()
|
|
c.draw_text(crs.to_proj4(), (50,10))
|
|
|
|
# save
|
|
c.save("testrenders/"+savename+".png")
|
|
|
|
|
|
|
|
|
|
|
|
|
|
# Source string generator
|
|
def sourcestrings(format):
|
|
# commonly used projections on global scale
|
|
# from http://www.remotesensing.org/geotiff/proj_list/
|
|
|
|
##Albers Equal-Area Conic
|
|
yield pycrs.utils.crscode_to_string("sr-org", 62, format)
|
|
##Azimuthal Equidistant
|
|
yield pycrs.utils.crscode_to_string("esri", 54032, format)
|
|
##Cassini-Soldner
|
|
# ...ignore, too specific
|
|
##Cylindrical Equal Area
|
|
yield pycrs.utils.crscode_to_string("sr-org", 8287, format)
|
|
##Eckert IV
|
|
yield pycrs.utils.crscode_to_string("esri", 54012, format)
|
|
##Eckert VI
|
|
yield pycrs.utils.crscode_to_string("esri", 54010, format)
|
|
##Equidistant Conic
|
|
yield pycrs.utils.crscode_to_string("esri", 54027, format)
|
|
##Equidistant Cylindrical
|
|
yield pycrs.utils.crscode_to_string("epsg", 3786, format)
|
|
##Equirectangular
|
|
yield pycrs.utils.crscode_to_string("sr-org", 8270, format)
|
|
##Gauss-Kruger
|
|
# ...not found???
|
|
##Gall Stereographic
|
|
yield pycrs.utils.crscode_to_string("esri", 54016, format)
|
|
##GEOS - Geostationary Satellite View
|
|
yield pycrs.utils.crscode_to_string("sr-org", 81, format)
|
|
##Gnomonic
|
|
# ...not found
|
|
##Hotine Oblique Mercator
|
|
yield pycrs.utils.crscode_to_string("esri", 54025, format)
|
|
##Krovak
|
|
yield pycrs.utils.crscode_to_string("sr-org", 6688, format)
|
|
##Laborde Oblique Mercator
|
|
# ...not found # yield pycrs.utils.crscode_to_string("epsg", 9813, format)
|
|
##Lambert Azimuthal Equal Area
|
|
yield pycrs.utils.crscode_to_string("sr-org", 28, format)
|
|
##Lambert Conic Conformal (1SP)
|
|
# ...not found
|
|
##Lambert Conic Conformal (2SP)
|
|
yield pycrs.utils.crscode_to_string("sr-org", 29, format) # yield pycrs.utils.crscode_to_string("epsg", 9802, format)
|
|
##Lambert Conic Conformal (2SP Belgium)
|
|
# ...ignore, too specific
|
|
##Lambert Cylindrical Equal Area
|
|
yield pycrs.utils.crscode_to_string("sr-org", 8287, format)
|
|
##Mercator (1SP)
|
|
yield pycrs.utils.crscode_to_string("sr-org", 16, format)
|
|
##Mercator (2SP)
|
|
yield pycrs.utils.crscode_to_string("sr-org", 7094, format)
|
|
##Miller Cylindrical
|
|
yield pycrs.utils.crscode_to_string("esri", 54003, format)
|
|
##Mollweide
|
|
yield pycrs.utils.crscode_to_string("esri", 54009, format)
|
|
##New Zealand Map Grid
|
|
# ...ignore, too specific
|
|
##Oblique Mercator
|
|
yield pycrs.utils.crscode_to_string("esri", 54025, format)
|
|
##Oblique Stereographic
|
|
yield pycrs.utils.crscode_to_string("epsg", 3844, format)
|
|
##Orthographic
|
|
yield pycrs.utils.crscode_to_string("sr-org", 6980, format)
|
|
##Polar Stereographic
|
|
yield pycrs.utils.crscode_to_string("sr-org", 8243, format)
|
|
##Polyconic
|
|
yield pycrs.utils.crscode_to_string("esri", 54021, format)
|
|
##Robinson
|
|
yield pycrs.utils.crscode_to_string("esri", 54030, format)
|
|
##Rosenmund Oblique Mercator
|
|
# ...not found
|
|
##Sinusoidal
|
|
yield pycrs.utils.crscode_to_string("sr-org", 6965, format)
|
|
##Swiss Oblique Cylindrical
|
|
# ...ignore, too specific
|
|
##Swiss Oblique Mercator
|
|
# ...ignore, too specific
|
|
##Stereographic
|
|
yield pycrs.utils.crscode_to_string("sr-org", 6711, format)
|
|
##Transverse Mercator
|
|
# ...not found???
|
|
##Transverse Mercator (Modified Alaska)
|
|
# ...ignore, too specific
|
|
##Transverse Mercator (South Oriented)
|
|
# ...ignore, too specific
|
|
##Tunisia Mining Grid
|
|
# ...ignore, too specific
|
|
##VanDerGrinten
|
|
yield pycrs.utils.crscode_to_string("sr-org", 6978, format)
|
|
|
|
# bunch of randoms
|
|
#yield pycrs.utils.crscode_to_string("esri", 54030, format)
|
|
#yield pycrs.utils.crscode_to_string("sr-org", 7898, format)
|
|
#yield pycrs.utils.crscode_to_string("sr-org", 6978, format)
|
|
#yield pycrs.utils.crscode_to_string("epsg", 4324, format)
|
|
#yield pycrs.utils.crscode_to_string("sr-org", 6618, format)
|
|
#yield pycrs.utils.crscode_to_string("sr-org", 22, format)
|
|
#yield pycrs.utils.crscode_to_string("esri", 54031, format)
|
|
# add more...
|
|
|
|
# Misc other crs for testing
|
|
#crs = pycrs.utils.crscode_to_string("esri", 54030, "proj4")
|
|
#crs = pycrs.utils.crscode_to_string("sr-org", 6978, "proj4")
|
|
#crs = pycrs.parser.from_sr_code(7898)
|
|
#crs = pycrs.parser.from_epsg_code(4324)
|
|
#crs = pycrs.parser.from_sr_code(6618)
|
|
#crs = pycrs.parser.from_sr_code(22)
|
|
#crs = pycrs.parser.from_esri_code(54031)
|
|
#proj4 = "+proj=longlat +ellps=WGS84 +datum=WGS84"
|
|
#proj4 = "+proj=aea +lat_1=24 +lat_2=31.5 +lat_0=24 +lon_0=-84 +x_0=400000 +y_0=0 +ellps=GRS80 +units=m +no_defs "
|
|
#proj4 = "+proj=larr +datum=WGS84 +lon_0=0 +lat_ts=45 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
|
|
#proj4 = "+proj=nsper +datum=WGS84 +ellps=WGS84 +lon_0=-60 +lat_0=40 +h=2000000000000000000000000"
|
|
|
|
|
|
|
|
|
|
|
|
# Testing format outputs
|
|
def testoutputs(crs):
|
|
print("ogc_wkt:\n")
|
|
try:
|
|
print(crs.to_ogc_wkt()+"\n")
|
|
global ogcwkt_outputs
|
|
ogcwkt_outputs += 1
|
|
except: logging.warn(traceback.format_exc())
|
|
|
|
print("esri_wkt:\n")
|
|
try:
|
|
print(crs.to_esri_wkt()+"\n")
|
|
global esriwkt_outputs
|
|
esriwkt_outputs += 1
|
|
except: logging.warn(traceback.format_exc())
|
|
|
|
print("proj4:\n")
|
|
try:
|
|
print(crs.to_proj4()+"\n")
|
|
global proj4_outputs
|
|
proj4_outputs += 1
|
|
except: logging.warn(traceback.format_exc())
|
|
|
|
|
|
|
|
|
|
|
|
#############################################################################
|
|
#############################################################################
|
|
#############################################################################
|
|
|
|
|
|
|
|
|
|
###########################
|
|
# From OGC WKT
|
|
print("--------")
|
|
print("Testing from ogc wkt:")
|
|
print("")
|
|
|
|
totals = 0
|
|
loaded = 0
|
|
ogcwkt_outputs = 0
|
|
esriwkt_outputs = 0
|
|
proj4_outputs = 0
|
|
renders = 0
|
|
|
|
for wkt in sourcestrings("ogcwkt"):
|
|
totals += 1
|
|
|
|
# test parsing
|
|
try:
|
|
print("From:\n")
|
|
print(wkt)
|
|
print("")
|
|
crs = pycrs.parser.from_ogc_wkt(wkt)
|
|
loaded += 1
|
|
|
|
# test outputs
|
|
print("To:\n")
|
|
testoutputs(crs)
|
|
|
|
# test render
|
|
try:
|
|
print("Rendering...")
|
|
savename = "%i_from_ogcwkt" % totals
|
|
render_world(crs, savename)
|
|
renders += 1
|
|
print("Successully rendered! \n")
|
|
except:
|
|
logging.warn(traceback.format_exc()+"\n")
|
|
|
|
except:
|
|
logging.warn(traceback.format_exc()+"\n")
|
|
|
|
print("Summary results:")
|
|
print(" Loaded: %f%%" % (loaded/float(totals)*100) )
|
|
print(" Outputs (OGC WKT): %f%%" % (ogcwkt_outputs/float(totals)*100) )
|
|
print(" Outputs (ESRI WKT): %f%%" % (esriwkt_outputs/float(totals)*100) )
|
|
print(" Outputs (Proj4): %f%%" % (proj4_outputs/float(totals)*100) )
|
|
print(" Renders: %f%%" % (renders/float(totals)*100) )
|
|
|
|
|
|
|
|
|
|
###########################
|
|
# From PROJ4
|
|
print("--------")
|
|
print("Testing from proj4:")
|
|
print("")
|
|
|
|
totals = 0
|
|
loaded = 0
|
|
ogcwkt_outputs = 0
|
|
esriwkt_outputs = 0
|
|
proj4_outputs = 0
|
|
renders = 0
|
|
|
|
for proj4 in sourcestrings("proj4"):
|
|
totals += 1
|
|
|
|
# test parsing
|
|
try:
|
|
print("From:\n")
|
|
print(proj4)
|
|
print("")
|
|
crs = pycrs.parser.from_proj4(proj4)
|
|
loaded += 1
|
|
|
|
# test outputs
|
|
print("To:\n")
|
|
testoutputs(crs)
|
|
|
|
# test render
|
|
try:
|
|
print("Rendering...")
|
|
savename = "%i_from_proj4" % totals
|
|
render_world(crs, savename)
|
|
renders += 1
|
|
print("Successully rendered! \n")
|
|
except:
|
|
logging.warn(traceback.format_exc()+"\n")
|
|
|
|
except:
|
|
logging.warn(traceback.format_exc()+"\n")
|
|
|
|
print("Summary results:")
|
|
print(" Loaded: %f%%" % (loaded/float(totals)*100) )
|
|
print(" Outputs (OGC WKT): %f%%" % (ogcwkt_outputs/float(totals)*100) )
|
|
print(" Outputs (ESRI WKT): %f%%" % (esriwkt_outputs/float(totals)*100) )
|
|
print(" Outputs (Proj4): %f%%" % (proj4_outputs/float(totals)*100) )
|
|
print(" Renders: %f%%" % (renders/float(totals)*100) )
|
|
|
|
|
|
|
|
|
|
###########################
|
|
# From ESRI WKT/PRJ FILE
|
|
print("--------")
|
|
print("Testing from esri wkt:")
|
|
print("")
|
|
|
|
totals = 0
|
|
loaded = 0
|
|
ogcwkt_outputs = 0
|
|
esriwkt_outputs = 0
|
|
proj4_outputs = 0
|
|
renders = 0
|
|
|
|
for wkt in sourcestrings("esriwkt"):
|
|
totals += 1
|
|
|
|
# test parsing
|
|
try:
|
|
print("From:\n")
|
|
print(wkt)
|
|
print("")
|
|
crs = pycrs.parser.from_esri_wkt(wkt)
|
|
loaded += 1
|
|
|
|
# test outputs
|
|
print("To:\n")
|
|
testoutputs(crs)
|
|
|
|
# test render
|
|
try:
|
|
print("Rendering...")
|
|
savename = "%i_from_esriwkt" % totals
|
|
render_world(crs, savename)
|
|
renders += 1
|
|
print("Successully rendered! \n")
|
|
except:
|
|
logging.warn(traceback.format_exc()+"\n")
|
|
|
|
except:
|
|
logging.warn(traceback.format_exc()+"\n")
|
|
|
|
print("Summary results:")
|
|
print(" Loaded: %f%%" % (loaded/float(totals)*100) )
|
|
print(" Outputs (OGC WKT): %f%%" % (ogcwkt_outputs/float(totals)*100) )
|
|
print(" Outputs (ESRI WKT): %f%%" % (esriwkt_outputs/float(totals)*100) )
|
|
print(" Outputs (Proj4): %f%%" % (proj4_outputs/float(totals)*100) )
|
|
print(" Renders: %f%%" % (renders/float(totals)*100) )
|
|
|
|
|
|
|