0

I programmed a fishnet in python to use as a reference map for georeferencing. I am facing the following issue:

The fishnet is not displayed correctly in ArcMap 10.5.1. The fishnet has parameters lat-values: 45 - 56, lon-values: 6 - 24 (in Ferro). The projection I use is custom made for the map I georeference it to which has a central meridian at 30° ferro (14° 20' in Greenwhich). Instead of displaying the map appropriately the coordinates of the fishnet are displayed at longitude of origin (30°) and latitude of origion (0°) (see Picture)enter image description here.

The fishnet is displayed correctly (in the sense that it has the correct coordinates) in QGIS. My best guess right now, is that it is a problem with the projection I am using. However, I am not sure what. Things I ruled out: it also does not work if I use the standard greenwhich as Pm.

The Code for the Fishnet:

    """
    Creating fishnet with lat-values: 45 - 56, lon-values: 6 - 24, graticule 1 x 1.
    """

   import numpy as np
   import geopandas as gpd
   from shapely.geometry import Polygon


   #Graticule of Fishnet
   sq_lat_deg = 0.5
   sq_lon_deg = 1 

   #Size of the fishnet and coordinates
   lats = np.arange(45, 56, sq_lat_deg)
   lons = np.arange(23, 42, sq_lon_deg)

   # Creating list of lat-lon values 
   coord_fishnet = []

   for lat in lats:
           for lon in lons:
                coord_fishnet.append(Polygon([(lon, lat), (lon, lat - 
                                              sq_lat_deg), (lon -
                                              sq_lon_deg, lat - 
                                              sq_lat_deg), (lon -
                                              sq_lon_deg, lat)]))

    # Setting Projection: Lambert conic projection, parallels 50/53, ellps bessel
    # datum? lat_0? Could set pm (prime meridian) to ferro.    
    crs = {"proj": "lcc", "lat_0": "0", "lon_0": "30", "lat_1": "50",
           "lat_2": "53", "ellps": "bessel", "datum": "potsdam", "pm": "ferro",
           "no_defs": True}

    # Creating ESRI readable shapefile
    fishnet = gpd.GeoSeries(coord_fishnet, crs =crs)

    fishnet.to_file(driver = 'ESRI Shapefile', filename = "fishnet.shp")

    # Testing
    fishnet.plot(color='white', edgecolor='black')

    test = gpd.read_file('fishnet.shp')
    test.plot(color='white', edgecolor='black') 

Your Answer

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

Browse other questions tagged or ask your own question.