Ticket #231: project_optimal_albers.py

File project_optimal_albers.py, 2.9 KB (added by bbest, 12 years ago)

Python ArcGIS code

Line 
1import arcgisscripting, os, sys
2
3# P:\data\gs\gs020924 P:\data\gs\gs020924a P:\data\gs\gs020924a.prj 1000
4
5# input variables
6in_geodataset = sys.argv[1] # r"C:\bbest\code\mget\MGET\Branches\HabMod\NMML_Lab\data\survey99-00.mdb\gshhs_extent"
7out_geodataset = sys.argv[2]
8out_prj = sys.argv[3]
9out_cellsize = sys.argv[4]  # for raster projection, defaults to "#"
10# "C:\bbest\NMML_Lab\data\shore\ply_extent_gcs.shp" "C:\bbest\NMML_Lab\data\shore\ply_extent_alb.shp" "C:\bbest\NMML_Lab\data\ply_extent_alb.prj" "#"
11
12
13gp = arcgisscripting.create()
14
15# get spatial geographic projection reference from input dataset
16d = {} # dictionary of variables for passing into projection string
17desc = gp.Describe(in_geodataset)
18sr = desc.SpatialReference
19d['GCSName'] = sr.GCSName
20d['DatumName'] = sr.DatumName
21d['SpheroidName'] = sr.SpheroidName
22d['SemiMajorAxis'] = sr.SemiMajorAxis
23d['InverseFlattening'] = 1/sr.Flattening
24d['PrimeMeridianName'] = sr.PrimeMeridianName
25d['Longitude'] = sr.Longitude
26d['AngularUnitName'] = sr.AngularUnitName
27d['RadiansPerUnit'] = sr.RadiansPerUnit
28# output units
29d['UnitName'] = 'Meter'
30d['UnitNumber'] = 1
31
32# x and y
33x1, y1, x2, y2 = [float(v) for v in desc.Extent.split()]
34d['x_0'] = x1 + (x2 - x1)/2 # Central_Meridian
35d['y_0'] = y1 + (y2 - y1)/2 # Latitude_Of_Origin
36dy = (y2-y1) / 6
37d['y_1'] = y1 + dy          # Standard_Parallel_1
38d['y_2'] = y2 - dy          # Standard_Parallel_2
39
40# write projection string
41f = open(out_prj, 'w')
42proj_str = '''PROJCS["Custom_Albers_Equal_Area_Conic",
43                     GEOGCS["%(GCSName)s",
44                            DATUM["%(DatumName)s",
45                                  SPHEROID["%(SpheroidName)s",%(SemiMajorAxis)1.1f,%(InverseFlattening)1.9f]],
46                            PRIMEM["%(PrimeMeridianName)s",%(Longitude)1.1f],
47                            UNIT["%(AngularUnitName)s",%(RadiansPerUnit)1.16f]],
48                     PROJECTION["Albers"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",%(x_0)g],
49                     PARAMETER["Standard_Parallel_1",%(y_1)g],PARAMETER["Standard_Parallel_2",%(y_2)g],PARAMETER["Latitude_Of_Origin",%(y_0)g],
50                     UNIT["%(UnitName)s",%(UnitNumber)g]]''' % d
51# TO CHECK: consistent floating point decimal places? otherwise doesn't recognize as same datum and requires a datum transformation
52f.write(proj_str)
53f.close()
54
55# project
56try:
57    if desc.DataType == 'RasterDataset':
58        #gp.ProjectRaster_management(in_geodataset, out_geodataset, out_prj, 'NEAREST', out_cellsize, geographic_transform, Registration_Point, in_coor_system)
59        gp.ProjectRaster_management(in_geodataset, out_geodataset, out_prj, 'NEAREST', out_cellsize) # , '#', '#', sr) #, geographic_transform, Registration_Point, in_coor_system)   
60    else:  # presume desc.DataType == 'FeatureClass'
61        gp.Project(in_geodataset, out_geodataset, out_prj)
62except:
63    print gp.GetMessages(2)
64