| 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 |
|
|---|
