RiccardoDC commented 5 years ago

Hello, When I run Code1 I get a reprojection that is slightly shifted (Output1). When I run Code 2 the reprojection is correct (Output2). I am curious to know why; I changed the input of the to_crs function from manhole.crs to sewer.crs. However, they have the same crs; in fact, they are perfectly aligned in both cases. It gets even worse if I try to reproject the sewer and the manhole using the crs of the buildings. Could you please help? Thank you very much, Riccardo

Code1: import geopandas as gpd

import files from openstreetmap

bld = gpd.read_file("buildings.shp") # buildings

import from dxf conversion

sewer= gpd.read_file("prova-line.shp") # sewer

import from R with CRS

manhole = gpd.read_file("junctions.shp") # manholes

shows coordinates system

print(bld.crs) print(sewer.crs) print(manhole.crs)

reprojection

bld = bld.to_crs(sewer.crs)

bld = bld.to_crs(manhole.crs)

creates centroids for buildings

bld['centroid_column'] = bld.centroid bld_centroid = bld.set_geometry('centroid_column')

figures

sewer_plot = sewer.plot(markersize = 4, figsize=(15,15), color = "black"); manhole.plot(ax=sewer_plot, markersize = 10, figsize=(15,15), color = "red"); bld.plot(ax=sewer_plot, markersize = 4, figsize=(15,15), color = "grey");

Output 1: {'init': 'epsg:4326'} {'init': 'epsg:21781'} {'proj': 'somerc', 'lat_0': 46.95240555555556, 'lon_0': 7.439583333333333, 'k_0': 1, 'x_0': 600000, 'y_0': 200000, 'ellps': 'bessel', 'units': 'm', 'no_defs': True}

Code 2:

reprojection

bld = bld.to_crs(sewer.crs)

bld = bld.to_crs(manhole.crs)

Output 2: {'init': 'epsg:4326'} {'init': 'epsg:21781'} {'proj': 'somerc', 'lat_0': 46.95240555555556, 'lon_0': 7.439583333333333, 'k_0': 1, 'x_0': 600000, 'y_0': 200000, 'ellps': 'bessel', 'units': 'm', 'no_defs': True}

phildias commented 5 years ago

Is it possible to provide access to the buildings.shp, prova-line.shp and junctions.shp files?

RiccardoDC commented 5 years ago

Thank you for answering. Sure! I attached the files. files.zip

ljwolf commented 5 years ago

Reading the .prj files of the three you've sent and also in the output you state above, none are actually precisely the same projection... What's the intended projection for each?

RiccardoDC commented 5 years ago

Hello, thanks for answering. Each files comes from a diffrent source. bld from openstreetmap, sewer from a dxt file and converted it as a shp file and, manhole from an inp file and converted with the inp_to_files function in R. The crs of the file sewer and manhole are the same, in fact they are always aligned in both cases. What I don't understand is why geopandas reprojects the bld correctly if uses the sewer projection while, incorrectly if it uses the manhole projection. I found this post (https://gis.stackexchange.com/questions/297528/wrong-reprojecting-of-a-shapefile-with-to-crs-in-geopandas) where they had the same curiosity but with no clear answer.

ljwolf commented 5 years ago

The files for sewer and manhole are close and have the same meaning, but they're indeed not identical files:

junctions.prj

PROJCS["Hotine_Oblique_Mercator_Azimuth_Center",
    GEOGCS["GCS_Bessel 1841",
        DATUM["D_unknown",
            SPHEROID["bessel",6377397.155,299.1528128]
        ],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]
    ],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",46.95240555555556],
    PARAMETER["longitude_of_center",7.439583333333333],
    PARAMETER["azimuth",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",200000],
    UNIT["Meter",1]
]

prova-line.prj

PROJCS["CH1903_LV03",
    GEOGCS["GCS_CH1903",
        DATUM["D_CH1903",
            SPHEROID["Bessel_1841",6377397.155,299.1528128]
        ],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]
    ],
    PROJECTION["Hotine_Oblique_Mercator_Azimuth_Center"],
    PARAMETER["latitude_of_center",46.95240555555556],
    PARAMETER["longitude_of_center",7.439583333333333],
    PARAMETER["azimuth",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",600000],
    PARAMETER["false_northing",200000],
    UNIT["Meter",1]
]

buildings.prj

GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137,298.257223563]
    ],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]
]

This looks like it bears similarity with the question you link on stackexchange because of the datum differences, although those look superficial... will need to dig.

RiccardoDC commented 5 years ago

Thank you very much! It is more clear to me now:) Thanks a lot.

jorisvandenbossche commented 5 years ago

@RiccardoDC @ljwolf this can be closed? (didn't fully follow the issue)

ljwolf commented 5 years ago

I'd leave open because it's:

replicable
clearly providing a wrong result
uncertain whether it's a pyproj or input data issue.