Georeferencing

2792
11
09-12-2016 11:14 AM
OliverLeimer1
New Contributor

Hello community!

I have a special question on georeferencing, on which I haven't yet found an answer.

I want to create danger areas for groundwater-pumps. To achieve that, I have a point-shapefile of all the groundwater-pumps in my city, including information of their pumping capacity, the groundwater-flow direction in that area, etc. 


In addition to that I have created an image in a program called "Processing Shemat", showing the spatial expansion of the pump-related thermal influence in the ground water from one single pumping station. This image was created for one single pumping station, but will look the same for all others, too.  So let's assume, that this one graphic is valid for all the points in my shapefile.

What I now want to do is importing that image  to arcgis and not only connect it with one of the points in my shapefile, but all of the points (so that as result I'll have 1000 images (one per point in my shapefile) instead of only one)?

All the tutorials on georeferencing in the internet showed me, how I can adjust that picture to a single point, but doing that for 1000 pumping stations in the city is too much of an effort. 
So is there a possibility of (sorry for my amateurish language) kind of telling arcgis to connect the image from one point in the picture to all the points in my shapefile?

Does that even make sense with a raster, or should I first use ArcScan to vectorize it and make a similar operation later?

Thanks in advance for your help  

0 Kudos
11 Replies
OliverLeimer1
New Contributor

The points rotation is based on the direction of the groundwater-flow, which in this area is relatively constant going to SE - ESE, referring to government data. 
I'm still not sure, whether I shall make the effort and calculate (validate) it by myself (which would require more data and would cause more waiting-time - the mills of my city government grind slowly).

0 Kudos
XanderBakker
Esri Esteemed Contributor

I just gave it a shot using one of the links that Rebecca Strauch, GISP provided. The result that I obtained is shown below:

What I did was:

  • I georeferenced the image with the polygons using the coordinates displayed (X from 0 - 100 and Y from 0 - 200).
  • I didn't convert the image to polygons, since that would generate a lot of polygons (the image is not very clear), but I captured the polygons manually (see attached gdb)
  • The polygons don't have the 0,0 where it should be (around X=49 and Y=185), so in the code I translate the polygons first to the real 0,0 so that the rotation will yield the desired results
  • All vertices of the polygons are read, corrected for the real 0,0 (translated  using X-49 and Y-185), then rotated using the angle (about 58° counter clockwise based on your map) and translated to each input point (stored in a list)

See the code below:

def main():
    import arcpy

    fc_pol = r'C:\GeoNet\Georef\data.gdb\Thermal_influence_int_sel'
    fc_out = r'C:\GeoNet\Georef\data.gdb\test04'
    sr = arcpy.Describe(fc_pol).spatialReference

    # center for rotation
    pnt_ori = arcpy.Point(49.435, 185.515)
    angle = 58 # counter clockwise

    # list of points
    lst_crds = [[613904.848, 5346348.604],
                [613825.804, 5346381.016],
                [614007.705, 5346382.008],
                [613867.145, 5346301.971],
                [613960.080, 5346242.109]]

    # initial rotation and translation
    lst_features = []
    with arcpy.da.SearchCursor(fc_pol, ('SHAPE@')) as curs:
        for row in curs:
            pol_in = row[0]
            for crds in lst_crds:
                lst_parts = []
                for part in pol_in:
                    lst_vertices = []
                    for vertice in part:
                        if vertice is None:
                            lst_vertices.append(vertice)
                        else:
                            # translate first (correct for 0,0 of polygons)
                            pnt1 = arcpy.Point(vertice.X - pnt_ori.X, vertice.Y - pnt_ori.Y)
                            # rotate and translate
                            tpl_out = trans(pnt1.X, pnt1.Y, crds[0], crds[1], angle)
                            pnt_out = arcpy.Point(tpl_out[0], tpl_out[1])
                            lst_vertices.append(pnt_out)
                    lst_parts.append(lst_vertices)
                pol_out = arcpy.Polygon(arcpy.Array(lst_parts), sr)
                lst_features.append(pol_out)

    # write polygons to output fc
    arcpy.CopyFeatures_management(lst_features, fc_out)


def trans(px, py, tx, ty, angle):
    # fuente: http://gis.stackexchange.com/questions/67425/rotating-and-moving-polygons-to-different-angles-using-...
    '''
    Rotate a point px,py around origin 0,0
    Sense clockwise or anticlockwise
      (default clockwise)
    Transform to plot coords tx,ty
    Does not handle z or m values
    '''
    import math
    import numpy
    sense = 'anticlock'

    # angle in degrees +/- 360
    th = math.radians(angle) # theta (radians)
    if sense == 'anticlock':
        matrix = [[math.cos(th),-math.sin(th),tx],
                  [ math.sin(th), math.cos(th),ty],
                  [0,0,1]]
    elif sense == 'clock':
        matrix = [[ math.cos(th),math.sin(th),tx],
                  [-math.sin(th),math.cos(th),ty],
                  [0,0,1]]
    else:
        raise Exception, "sense error in DrawPlotR"

    a = numpy.array(matrix)
    b = numpy.array([px,py,1])
    c = numpy.dot(a,b)
    return (c[0],c[1])

if __name__ == '__main__':
    main()