@@ -0,0 +1,1632 @@
+def my_project(in_east, in_north, direction):%0A    '''%0A    Sample user-defined projection and inverse projection of (lon,lat) to (x,y) %0A      %0A      function %5Bout_east,out_north%5D = my_project(in_east,in_north,direction) %0A      %0A      DESCRIPTION:%0A         Define projections between geographical and Euclidean coordinates %0A      %0A      INPUT: %0A        in_east   = 1D vector containing longitude (forward) x (reverse)%0A        in_north  = 1D vector containing latitude  (forward) y (reverse)%0A        direction = %5B'forward' ;  'inverse'%5D%0A                %0A      OUTPUT:%0A        (lon,lat) or (x,y) depending on choice of forward or reverse projection%0A      %0A      EXAMPLE USAGE%0A         lon,lat = my_project(x,y,'reverse') %0A    '''%0A    from mpl_toolkits.basemap import pyproj%0A#    import mpl_toolkits.basemap.pyproj as pyproj%0A    #state_plane = pyproj.Proj(r'+proj=tmerc +datum=NAD83 +lon_0=-70d10 lat_0=42d50 k=.9999666666666667 x_0=900000 y_0=0 +to_meter=1')%0A    #state_plane = pyproj.Proj(r'+proj=tmerc +lat_0=42.83333333333334 +lon_0=-70.16666666666667 +k=0.9999666666666667 +x_0=900000 +y_0=0 +ellps=GRS80 +units=m +no_defs')%0A    state_plane = pyproj.Proj(r'+proj=tmerc +lat_0=42d50 +lon_0=-70d10 +k=0.9999666666666667 +x_0=900000 +y_0=0 +ellps=GRS80 +units=m +no_defs')%0A    %0A    wgs = pyproj.Proj(proj='latlong', datum='WGS84', ellps='WGS84')%0A    if direction=='forward':%0A        lon = in_east%0A        lat = in_north%0A        x,y = pyproj.transform(wgs, state_plane, lon, lat)%0A        return x, y%0A    else:%0A        x = in_east%0A        y = in_north%0A%0A        lon, lat = pyproj.transform(state_plane, wgs, x, y)%0A        return lon, lat%0A

add my_project.py