1 | import arcgisscripting, os, sys
|
---|
2 |
|
---|
3 | # P:\data\gs\gs020924 P:\data\gs\gs020924a P:\data\gs\gs020924a.prj 1000
|
---|
4 |
|
---|
5 | # input variables
|
---|
6 | in_geodataset = sys.argv[1] # r"C:\bbest\code\mget\MGET\Branches\HabMod\NMML_Lab\data\survey99-00.mdb\gshhs_extent"
|
---|
7 | out_geodataset = sys.argv[2]
|
---|
8 | out_prj = sys.argv[3]
|
---|
9 | out_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 |
|
---|
13 | gp = arcgisscripting.create()
|
---|
14 |
|
---|
15 | # get spatial geographic projection reference from input dataset
|
---|
16 | d = {} # dictionary of variables for passing into projection string
|
---|
17 | desc = gp.Describe(in_geodataset)
|
---|
18 | sr = desc.SpatialReference
|
---|
19 | d['GCSName'] = sr.GCSName
|
---|
20 | d['DatumName'] = sr.DatumName
|
---|
21 | d['SpheroidName'] = sr.SpheroidName
|
---|
22 | d['SemiMajorAxis'] = sr.SemiMajorAxis
|
---|
23 | d['InverseFlattening'] = 1/sr.Flattening
|
---|
24 | d['PrimeMeridianName'] = sr.PrimeMeridianName
|
---|
25 | d['Longitude'] = sr.Longitude
|
---|
26 | d['AngularUnitName'] = sr.AngularUnitName
|
---|
27 | d['RadiansPerUnit'] = sr.RadiansPerUnit
|
---|
28 | # output units
|
---|
29 | d['UnitName'] = 'Meter'
|
---|
30 | d['UnitNumber'] = 1
|
---|
31 |
|
---|
32 | # x and y
|
---|
33 | x1, y1, x2, y2 = [float(v) for v in desc.Extent.split()]
|
---|
34 | d['x_0'] = x1 + (x2 - x1)/2 # Central_Meridian
|
---|
35 | d['y_0'] = y1 + (y2 - y1)/2 # Latitude_Of_Origin
|
---|
36 | dy = (y2-y1) / 6
|
---|
37 | d['y_1'] = y1 + dy # Standard_Parallel_1
|
---|
38 | d['y_2'] = y2 - dy # Standard_Parallel_2
|
---|
39 |
|
---|
40 | # write projection string
|
---|
41 | f = open(out_prj, 'w')
|
---|
42 | proj_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
|
---|
52 | f.write(proj_str)
|
---|
53 | f.close()
|
---|
54 |
|
---|
55 | # project
|
---|
56 | try:
|
---|
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)
|
---|
62 | except:
|
---|
63 | print gp.GetMessages(2)
|
---|
64 |
|
---|