New answers tagged coordinate-system
0
I think I figured out the root of the issue with help from a friend--since Albers is a regional projection, its parameters aren't valid for a global extent. To get the raster into the Albers projection, I first clipped it down to an extent slightly larger than the land area that would need to be shown in the map and then reprojected it. That did the trick!
0
I solved the problem with memory driver and gdal Warp. I share my code and my resources in below.
def calculate_ndvi(self,red_path,nir_path,output_path):
g = gdal.Open(red_path)
red = g.ReadAsArray().astype('float')
g = gdal.Open(nir_path)
nir = g.ReadAsArray().astype('float')
numpy.seterr(divide='...
1
For a formula that you can implement, see this page for an explanation and scroll down the bottom for various implemenations 3857 to/from 4326
Python Implementation provided:
def MetersToLatLon( mx, my ):
"Converts XY point from Spherical Mercator EPSG:900913 to lat/lon in WGS84 Datum"
lon = (mx / originShift) * 180.0
lat = (my / originShift) ...
4
I see that you put a +towgs84 parameter in your destination crs (WGS84), but it should be included in the source crs or in the transformation (as dx, dy, dz parameters), not the destination crs.
If I transform from:
+proj=lcc +lat_1=44.33333333333334 +lat_2=46 +lat_0=43.66666666666666 +lon_0=-120.5 +x_0=609601.2192024384 +y_0=0 +ellps=clrk66 +towgs84=-8,...
1
I would recommend either:
pandas and pyproj.Transforner for the most efficient method:
import pandas
df = pandas.read_csv(...)
Then see: https://gis.stackexchange.com/a/334307/144357
Use geopandas with the .to_crs() method.
https://geopandas.readthedocs.io/en/latest/projections.html
import geopandas
import pandas
df = pandas.read_csv(...)
gdf = df....
0
Just a shot in the dark here as I am not able to reproduce the code you have provided at the moment.
The error suggests you should try to open the file in text mode. From the Python documentation on the open() function:
The first argument is a string containing the filename. The second
argument is another string containing a few characters describing ...
2
Are you shure that this point lies on that border? Intersecting points with lines is like spliting hair with an axe. It is possible, but it is very hart... You have to deal with 64 bit floating point precision coordinates and if there is marginal difference between them, there is no intersection. You should try to search polygons within small distance around ...
3
Using pyproj.CRS you can get the geodetic CRS of the projected CRS.
from pyproj import CRS
crs = CRS(infile.crs)
print(repr(crs.geodetic_crs))
<Geographic 2D CRS: GEOGCRS["GCS_Moon",DATUM["D_Moon",ELLIPSOID["Moon" ...>
Name: GCS_Moon
Axis Info [ellipsoidal]:
- lon[east]: Longitude (degree)
- lat[north]: Latitude (degree)
Area of Use:
- undefined
...
10
Does qgis just pick a default projection or does that mean that in epsg 4326 there's a 2d projection defined as well?
QGIS doesn't project the geographic coordinates to the map canvas. If the map has defined a geographic coordinate reference system, it gives to the X axis the angular coordinates of the Longitude axis, and to the Y axis the angular ...
2
For #1, yes, they are the same.
For #2, the UTM projection is a system of zones setup for the Transverse Mercator projection. See: https://en.m.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system
The string you are using to define the UTM projection is a PROJ string an each parameter is defined here: https://proj.org/operations/projections/...
4
The coordinates are based on the NAD27 datum, converted to UTM zone 11N.
You can use the TAGIS Converter to see them and get the WGS 84 geographic coordinates.
1
I am not sure if you considered that if your warp operation is rotating the image then the original top-left pixel is not at the top-left corner of the envelope of the warped image.
A test with gdalwarp and image LC08_L1TP_033033_20180706_20180717_01_T1_B5.tif
Command
gdalwarp -t_srs epsg:4326 LC08_L1TP_033033_20180706_20180717_01_T1_B5.tif warped.tif -co ...
3
Try printing out the "dsc.spatialReference.Name". When I did, none of my Spatial Reference Names had spaces and instead had underscores "_". I suspect that may be why none of your spatial references match "NAD 1983 StatePlane California III FIPS 0403 (US Feet)".
So it instead should be "NAD_1983_StatePlane_California_III_FIPS_0403_(US_Feet)".
0
A translation of coordinates is an affine conversion.
It can be done, with pyproj, through a PROJ Affine Transformation pipeline:
from pyproj import Transformer
import numpy as np
np.set_printoptions(precision=3,suppress=True)
# points list in [X,Y,Z]
points = [[0,0,0],
[1,0,0],
[0,1,0],
[0,0,1],
[-1,0,0],
[0,-1,...
0
I'm running QGis 3.10.2 on Ubuntu 18.04.4 LTS
QGIS version 3.10.2-A Coruña
QGIS code revision 616ad4531b
Compiled against Qt 5.9.5
Running against Qt 5.9.5
Compiled against GDAL/OGR 2.2.3
Running against GDAL/OGR 2.2.3
Compiled against GEOS 3.6.2-CAPI-1.10.2
Running against GEOS 3.7.1-CAPI-1.11.1 27a5e771
Compiled against SQLite 3.22.0
Running against ...
1
AAEXs, AAEYs and AAEZs are not any absolute or relative coordinates, but the (observed/recorded) directions of the magnetic field in the X, Y and Z directions (expressed) in nano Tesla (nT).
The measurement device (i.e., "sensor") was placed at 09.033 North 038.767 East and 2442 meter above the Mean Sea Level (most probably EGM2008).
1
This case of referencig is linear, because you are referencing a rectangular grid (row, pixel) to a rectangular grid (latitude, longitude).
The easiest way to do it, is with a World File.
Create a new text document and paste the following lines in it:
0.18
0.0
0.0
-0.18
-179.91
89.91
Save it as venus_shaded.wld, in the same folder as venus_shaded....
0
This smells like GDAL and proj installation/versioning issues.
QGIS uses GDAL/OGR to input and output many types of layers, and the proj library to handle projections.
QGIS version 3.10 transitioned to proj version 6.x. This has been a major change, with glitches appearing and being fixed in various point releases of QGIS, GDAL, and proj released in a ...
3
Not sure if transforming vector is worth it, but here's the trick.
According to docs, .transform() is ony applied to an ee.Feature, not ee.FeatureCollection. So you have to map the procedure over a FeatureCollection.
// Define an arbitrary region in which to compute random points.
var region = ee.Geometry.Rectangle(-119.224, 34.669, -99.536, 50.064);
// ...
0
The user running the query does not have public in its path.
The 1st query does not throw a not found error because the path is explicitly written (SELECT public.AddGeometryColumn ...)
You can check if the schema where PostGIS is installed is in the user search path by issuing
show search_path;
If not, you can permanently add the path by altering the ...
0
When using pyproj, note the differences from various releases in how it is used to transform data. Here are a few examples using new/old capabilities based on the question:
Using pyproj >= 2.2.0
import pyproj
print(pyproj.__version__) # 2.4.1
print(pyproj.proj_version_str) # 6.2.1
proj = pyproj.Transformer.from_crs(3857, 4326, always_xy=True)
x1, y1 = (-...
1
In R you can create a grid, please take a look at this documentation: https://www.rdocumentation.org/packages/mapplots/versions/1.5.1/topics/make.grid
Please take a look at this page, that might be what you need: (images taken form the same page, attributed to the author "Huanfa Chen")
https://rpubs.com/huanfaChen/grid_from_polygon
I hope this helps.
-1
As commented by @Vince:
Never, ever use the Polygon constructor without a SpatialReference
(spatial_reference is the second parameter). But your problem is
failure to create new Point objects, so the Array just has one
repeating vertex. The code samples provide the template, and they
create new Point objects for each loop iteration (it's best ...
3
Upon inspecting the dataset, I realized that the units of the data are in radians.
import xarray
import rioxarray
xds = xarray.open_dataset("OR_ABI-L2-LSTF-M6_G17_s20200341900321_e20200341909388_c20200341910038.nc")
Inside the x variable the attributes say the data is in radians:
xds.x.attrs
{'units': 'rad',
'axis': 'X',
'long_name': 'GOES Projection ...
1
The problem seems to be the definition of a Derived Engineering CRS with a Base Geodetic CRS.
As seen in the Section 14.6 (Derived engineering CRS) of the Geographic information — Well-known text representation of coordinate reference systems standarization (2.0.6 version), a Derived Engineering CRS must have a Base Engineering CRS, wich must have an ...
1
It is unusual to be trying to define the longitude from 0 to 359. Longitude is customarily referred to as either -180 to +180 or 180W to 180E.
That being said, you should investigate the Projections and Transformations Toolbox in ArcGIS. If you already have a layer in your desired output projection you can use Project tool. If you don't already have another ...
2
OSM nodes coordinates have 7 decimals. PostGIS stores geometries with 15 decimals. The difference between the two is noise and should not be relied on. In fact, you should instruct st_asText to show only 7 digits.
select st_astext(st_transform(way, 4326),7)
from planet_osm_line
where osm_id = 482283890;
There is an interesting note in the doc, hidden the ...
4
The tangent plane CRS is known as a Topocentric CRS.
The conversion formulae from Geocentric, and from Geographic coordinates to Topocentric coordinates are defined in the Section 4.1.2, of the IOGP document: Geomatics Guidance Note 7, part 2.
The EPSG registry has two example definitions of Topocentric CRS: EPSG:5819 and EPSG:5820.
For instance, we ...
2
The question was missing a little bit of info but with the dataset in hand, I was able to roughly find out what was going on.
First off, each point is actually composed of 6 different points that each contain a single attribute that belongs to the actual point, which is the one in the middle:
You need to isolate the true point with the tool Extract by ...
1
I found a Python Library called utm 0.5.0.
The conversion.py module in that library seems to do the trick.
I've modified the code by removing the reference to numpy and the custom OutOfRangeError library.
import math as mathlib
__all__ = ['to_latlon', 'from_latlon']
K0 = 0.9996
E = 0.00669438
E2 = E * E
E3 = E2 * E
E_P2 = E / (1.0 - E)
SQRT_E = ...
0
Judging from the world file, I assume your image represents the whole world (from -180° to 180° longitude and 90° to -90° latitude), then if you know your resized image dimensions (let's notate them imgsize_x,imgsize_y), you can calculate the long/lat position of a pixel center like this:
longitude = ((pixel_x - 0.5) * 360 / imgsize_x) - 180
latitude = 90 - ...
6
The main problem is in the image warping parameters, they are not giving you the right gridded image.
About the projection itself, you are using the Lambert Conical projection in the right way. You have just not information about the datum transformation.
A datum transformation is going to be necessary if we want to compare this map with one based on ...
4
Envelopes are often very inaccurate, depending on the coordinate reference system (CRS), or converting an envelope between two CRS. An envelope is 2 points, so projecting or transforming them to a different CRS means you can miss distortions occurring at the other two corners, on the edges, etc. An envelope in one CRS, if it was densified first, many times, ...
0
Have you tried the Define Projection tool in ArcGIS? It doesn't modify anything to the original file, just updates the projection information!
2
Looking at this page which proposes this proj string for this Italy South grid:
+proj=lcc +lat_0=39.5 +lat_1=39.5 +lon_0=14 +k_0=0.99906 +x_0=700000
+y_0=600000 +ellps=bessel +units=m +no_defs
Seems to indicate you are on the right track, that may be as accurate as the projection gets?
1
Except with OpenLayers, I don't know about out of the box solutions to consume Spherical Mercator tiles and display them with another projection.
You should use Mapproxy https://mapproxy.org (a server side solution) to consume the tiles and serve them while reprojecting them. You can make WMS calls and behind the scene, the tool will stitch the tiles
0
As @Vince suggested, Google Maps uses WGS84.
And as @mkennedy pointed out, the SRID of WGS84 is 4326 (epsg-registry.org).
Therefore, the SQL to generate the hyperlink/URL would be:
select
'http://maps.google.com/maps?q=' || chr(38) || 'layer=c' || chr(38) || 'cbll='
|| latitudey_dd || ',' || longitudex_dd as street_view_url
from
(
select
...
4
Problem with “On the Fly Projection” (OTF) in QGIS 3.4 when upgrading from QGIS version 2.8 [closed]
About the QGIS 2 On-The-Fly reprojection, in QGIS 3:
In QGIS 3, in Project ~> Properties ~> CRS, you can define the Coordinate Reference System (CRS) for the project map canvas.
If the canvas has a defined CRS, then QGIS will reproject the layers on-the-fly to that CRS before show them on the map.
If "No projection" is ticked on, then the project map ...
2
As a quick answer, most transforms are more complex than a scale/translate/unit, so you are looking for exceptions. The distortion come from the change in datum (that is the transformation, from lat/long to lat/long) and from the type of projection (from lat/long to XY).
Having the same (or rotated around the vertical) datum is a prerequisite but, if you ...
1
I don't really get why you can't use Geography data directly, and let Postgis do the work. This kind of computation usually work in a few ms, the computation itself is generally way less than or the other costs (like reading data). If you have performance issues, maybe the problem is not the projection but the request or the indexes.
If you really are in a ...
0
Regarding the use of HDF format with GDAL the key concept is that the file contains subdatasets as opposed to raster bands.
With a GeoTiff, for example, you can open it, and pick a band to read
f = ''
ds = gdal.Open(f)
if ds is None:
raise Exception('Failed to open dataset')
b = ds.GetRasterBand(1)
# etc...
When using a HDF file, the file itself ...
2
Proj6 seems to consider that EPSG:4258 and EPSG:4326 mean the same and conversion between them does not change anything.
projinfo -o PROJ -s "EPSG:4258" -t "EPSG:4326"
Candidate operations found: 1
-------------------------------------
Operation n 1:
EPSG:1149, ETRS89 to WGS 84 (1), 1.0 m, Europe - ETRS89
projinfo -k operation EPSG:1149
PROJ string:
+proj=...
6
I tried several coordinate systems for Florida and the various State Plane Florida East zone but the coordinates just didn't seem to fit.
Google kept finding the address in Hempstead, NY, so I tried the Long Island zone, unit of US survey feet, and the results look good. Try EPSG::6539 for the NAD83 (2011) version or EPSG::2908 for NAD83 (HARN). For the ...
1
You should be able to solve this by creating a src/main/resources/org/geotools/referencing/factory/gridshift/ folder (it will get rolled up into your jar when you do a mvn install) to drop the GIS_GRID.gsb file into.
You can use some code like this to check the transform is being picked up:
CoordinateReferenceSystem epsg31255 = CRS.decode("epsg:31255");
...
0
Some additional information:
To get the code in Ian Turton's answer (GeoTools 22.2, Java 8) to run properly I had to add the following dependencies to the pom.xml file:
<dependency>
<groupId>org.geotools</groupId>
<artifactId>gt-geotiff</artifactId>
<version>${geotools.version}</version>
</dependency&...
4
The +init=epsg:xxxx syntax (or in form of a {'init': 'epsg:xxx'}) is being deprecated in PROJ / pyproj (see here. You are indeed not using it directly, but, the previous releases of GeoPandas (<= 0.6) were still generating this format under the hood (e.g. when reading files, as you can see from the output of df.crs). That is the reason you see this ...
1
Why is this happening?
Because lot of changes in the way that QGIS handle datum transformations are being made, and seems like in your version (I have the same) some one of those changes has introduced a bug.
This changes are necessary because lot of changes in PROJ are being made.
Take into account that PROJ.4 strings are not the best way to ...
0
Your WKT doesn't define any axis order, just lists them.
You could try with WKT 2 that defines axes order like:
GEODCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1.0]]],
CS[ellipsoidal,2],
AXIS["latitude",north,ORDER[1]],
AXIS["longitude",east,ORDER[2]],
ANGLEUNIT["degree",...
0
Let's assume that the Earth is flat and that the meridians are parallel. Let's assume also that BHL-No.3080 point has coordinates Easting = 0 , Northing = 0. We need it to calculate the coordinates of the other points in a planimetric way.
Here we go: We will make a spreadsheet for the calculation of the planimetric coordinates:
The first row has the ...
3
If you have a matching .wld and .prj file then the gt-image module can read and write jpeg and png files.
Once you've read the file in to a GridCoverage you can reproject it and write it out with no difficulty:
AbstractGridFormat format = GridFormatFinder.findFormat(input);
Hints hints = null;
if (format instanceof GeoTiffFormat) {
hints = new Hints(...
Top 50 recent answers are included
