2

I have two GeoPandas GeoDataFrame objects with their corresponding geometry columns.

  1. GeoDataFrame 1: counties --> has polygons and multipolygons of Dominican Republic counties. The source data was obtained from the Dominican Republic 2010 census website and consist of shape files. It was apparently made using ArcGis.
  2. GeoDataFrame 2: hospitals --> contains longitude and latitude coordinates (Points) for medical centers in the Dominican Republic. The source data was obtained from the Dominican Health Ministry website and consist of a Microsoft Excel file. The coordinates appears to align well with Google Maps.

I am trying to overlay the hospitals (Points) on top of the counties (Polygons). The problem I have is that the two geometries don't match. The hospitals are in longitude/latitude coordinates, but I'm not sure about the counties coordinates.

How can I make the two coordinates (X and Y axis) match so that I can overlay them?

The easiest way to understand my problem is to take a look at the plots:

https://drive.google.com/open?id=14xYmNJh958q1XANdq7dz98wn1Nrxymyi

Some sample code below. Note the total_bounds attribute.

import pandas as pd
import geopandas as gpd

# Import counties data
counties = gpd.read_file('ShapeFilesCenso2010/MUNCenso2010.shp')
counties.crs = {'init':'epsg:3395'}

# Import hospital coordinates
hospitals = pd.read_excel('hospitals.xlsx')

# After processing the excel file I turn it into the geometry below:
print(hospitals.geometry.head())

0    POINT (-69.89062 18.49381)
1    POINT (-69.88906 18.47107)
2    POINT (-69.89158 18.49421)
3    POINT (-69.92358 18.45247)
4    POINT (-69.86830 18.50938)
Name: geometry, dtype: geometry

print(hospitals.crs)

{'init': 'epsg:3395'}

print(hospitals.geometry.total_bounds)

[-72.      ,  17.      , -68.4513  ,  19.884642]

# The output from the counties geometry
print(counties.geometry.head())

0    POLYGON ((397122.651 2050791.921, 397142.013 2...
1    MULTIPOLYGON (((298172.939 2034259.341, 298374...
2    POLYGON ((333950.281 2046005.250, 334156.656 2...
3    MULTIPOLYGON (((300647.656 2050014.250, 300638...
4    POLYGON ((303931.923 2102011.319, 304050.231 2...
Name: geometry, dtype: geometry

print(counties.crs)

{'init': 'epsg:3395'}

print(counties.geometry.total_bounds)

[182215.7656, 1933511.9638,  571429.3273, 2205216.25]

ADDENDUM:

Checking an XML file that came with the counties (census) data, I see the following attribute and wonder if it has anything to do with the problem:

<Process Name="DefineProjection_2" ToolSource="C:\Program Files\ArcGIS\ArcToolbox\Toolboxes\Data Management Tools.tbx\DefineProjection" Date="20080220" Time="081209">
DefineProjection D:\DIVTER_ENERO2008\dm_1207_v1.shp PROJCS['WGS_1984_UTM_Zone_19N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['False_Easting',500000.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-69.0],PARAMETER['Scale_Factor',0.9996],PARAMETER['Latitude_Of_Origin',0.0],UNIT['Meter',1.0]] D:\DIVTER_ENERO2008\dm_1207_v1.shp
</Process>
  • Based on your comment, I defined the hospital points with epsg:4326. I then used hospitals.to_crs(epsg=3395). This seems to be much better because the Y axis of both GeoDataFrames are roughly in the same range. However, the X axis if off. please see this plot to get what I mean: new plot – Carlos Muñiz Dec 15 '19 at 4:17
  • Checking an XML file that came with the counties (census) data, I see the following attribute and wonder if it has anything to do with the problem (please my edit to the original post above because the XML attribute is too long to insert here). – Carlos Muñiz Dec 15 '19 at 4:20
  • 1
    It's good you put the extra information in an edit. Comments are for requesting more info, not for adding it. – user2856 Dec 15 '19 at 4:35
2

The CRS of your hospital points is incorrect. EPSG:3395 (WGS 84 / World Mercator) uses meters for coordinates. Redefine it as EPSG:4326 (WGS84 / Geographic) which uses degrees (lon/lat) for coordinates.

The CRS of your counties polygons may be incorrect. The metadata you added shows the CRS is WGS 84 / UTM zone 19N (EPSG:32619) not EPSG:3395.

Once you fix the CRSs, you can reproject one of the layers to match the other and they should overlay.

counties.crs = {'init':'epsg:32619'}                  # set correct projection
hospitals.crs = {'init' :'epsg:4326'}                 # set correct projection
hospitals = hospitals.to_crs({'init': 'epsg:32619'})  # reproject to match

You may not even need to set the CRS for the counties layer, it's probably only wrong because you set it to the incorrect CRS (EPSG:3395) after reading it in:

# Import counties data
counties = gpd.read_file('ShapeFilesCenso2010/MUNCenso2010.shp')
counties.crs = {'init':'epsg:3395'}  # <=== Don't do this...
  • This worked excellent. I still have some minor modifications to do because it seems the Ministry of Health provided some wrong hospital coordinates. Also there may be a bit of a scaling mismatch. But this is the result so far: final result. Much appreciated! – Carlos Muñiz Dec 15 '19 at 5:35
  • @Carlos what do you mean "scaling mismatch"? – user2856 Dec 16 '19 at 0:32

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.