2

I was wondering if the following approach a) would be considered valid and b) someone already did this in python.

I have a python script where users can choose an area of analysis (anywhere on the world). The tool will select data based on this rectangle (data in WGS1984 lat/lon) and process it. The result is a shapefile (from Shapely) that can be mapped in QGIS or Arcgis (etc.). My current problem is: the users of my script don't have any background in GIS. Therefore, the tool should do as much automatically as is possible. As a further automation step, I would like to output the resulting Shapefile in a UTM-Projected Coordinate System suitable for the chosen area.

  • users will always choose a pretty small area (local scale), therefore, it is very unlikely that an area stretches significantly across two UTM-Zones
  • Minimal Distance/Area distortions are important, but only to some degree - the users will perform some calculations that require consistent distances (e.g. GetisOrd Statistic) UTM-Grid UTM-Grid on Wikimedia I think I should be able to select the best UTM-Zone by intersecting with the UTM-Grid:

Is there any python package that will detect the best UTM-Zone given an lat/lng point (and perhaps an area extent)? Are there any concerns regarding my approach?

3

There's the utm package, a bidirectional UTM-WGS84 converter for python. It's able to convert a (latitude, longitude) tuple into an UTM coordinate:

import utm
utm.from_latlon(51.2, 7.5)
>>> (395201.3103811303, 5673135.241182375, 32, 'U')
  • Many thanks Antonio! So I can get the UTM Grid Coordinate "32U" based on a lat/lng pair. How would I select the matching projected Coordinate System (e.g. by getting the SRID))? For example using: project = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='epsg:26913')) # destination coordinate system – Alex Jan 29 '18 at 9:59
  • I think I need the SRID to assign a projection system to my shapefile (example: Fiona): with fiona.open(shapeout, 'w', crs=from_epsg(3996), driver='ESRI Shapefile', schema=myschema) as output: [...] – Alex Jan 29 '18 at 10:05
  • Found the answer: Here's a (very similar) question and an answer that provides code for extracting the matching EPSG Code (SRID) based on UTM-Coordinates from Lat/Lng – Alex Jan 29 '18 at 10:18
1

Alright, the answer from Antonio above is definitely right and pointed me in the correct direction. Here is the complete code:

#convert_wgs_to_utm function, see https://stackoverflow.com/questions/40132542/get-a-cartesian-projection-accurate-around-a-lat-lng-pair
def convert_wgs_to_utm(lon, lat):
    utm_band = str((math.floor((lon + 180) / 6 ) % 60) + 1)
    if len(utm_band) == 1:
        utm_band = '0'+utm_band
    if lat >= 0:
        epsg_code = '326' + utm_band
    else:
        epsg_code = '327' + utm_band
return epsg_code
#Get Bounds of user selected area
bound_points_shapely = geometry.MultiPoint([(limLngMin, limLatMin), (limLngMax, limLatMax)])
get lat/lng from center
input_lon_center = bound_points_shapely.centroid.coords[0][0] #True centroid (coords may be multipoint)
input_lat_center = bound_points_shapely.centroid.coords[0][1]
#Get best UTM Zone SRID/EPSG Code for center coordinate pair
utm_code = convert_wgs_to_utm(input_lon_center, input_lat_center)
#define project_from_to with iteration
#see https://gis.stackexchange.com/questions/127427/transforming-shapely-polygon-and-multipolygon-objects
project = lambda x, y: pyproj.transform(pyproj.Proj(init='epsg:4326'), pyproj.Proj(init='epsg:{0}'.format(utm_code)), x, y)
#define schema
schema = {
    'geometry': 'Polygon',
}
with fiona.open('projectedShapefile.shp', mode='w', encoding='utf-8', driver='ESRI Shapefile', schema=schema,crs=from_epsg(utm_code)) as c:
    #for each shapely geometry:
    geom_proj = transform(project, geom)
    c.write({
    'geometry': geometry.mapping(geom_proj)})

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.