texts: ArcGIS tools say completed immediately but give no output

texts:
I have a problem working with ArcGIS. Sometimes (not always) when I run a geprocessing tool, I immediately get the pop-up with "completed" and a green tickmark. However, there is no output created, and it not realistic that the processing should finish so quickly.

I have not been able to discover a system as to when the tools work and when they do not. I have also tried with shorter names and putting them in locations with no spaces in the path string.

An example of a tool that does not work is Resample, but I also have problems with other tools.

Any suggestions?


texts: IFeatureIndex2.Index Performance Issues ArcObjects

texts:
I have a large number of polygon features in a personal geodatabase featureclass. I would like to find the nearest polygon for several points inside an area of interest(the red area in figure).

I found IIndexQuery2.NearestFeature is a good option since it 'Finds the nearest feature in index to the input shape' as per the documentation.

At present, I do not see any way of utilizing the existing spatial index in the geodatabase through arcobjects. The only option i found is to create a temporary index using IFeatureIndex2.Index method.

My problems are with the IFeatureIndex2.Index method. I only need to index a particular group of features. There is a IFeatureIndex2.FeatureCursor option  which supposedly reduces the number of features indexed. So i passed a featureCursor of the specified area to IFeatureIndex2.But mysteriously, this option only increases the time taken for index creation.

Also the Envelope parameter in the Index method (shown below) does not seem to have an effect on the performance.

public void Index (ITrackCancel pTrackCancel, IEnvelope pQueryingGeometryFullExtent));

It does not matter if null is passed or a valid envelope of Area is passed to this method as parameter. It takes same time to complete.

I tested it on a polygon featureclass with around 100000 features.
The Index Method took 4 seconds
Index Method with FeatureCursor took 12 seconds.

Can anyone suggest whether I am using the methods in the proper way? I might have to run this on an SDE with millions of features. So I am looking forward to indexing only the particular area so that the nearest feature queries would execute faster.


texts: Integer rasters referring to same scale using ArcMap

texts:
I'm using ArcMap 10. I want to display the area of simulated avalanches. I want this area to be divided in different colors depending on the thickness of the deposit. I managed to do that. But I want all the simulated avalanches to refer to the same scale and same legend of thickness. This means that for every simulation a certain range of values of thickness (e.g. from 2m to 4m) falls within the same color. The problem is I need to transform my raster into polygons (for future calculations of area ratios) and to do this it must be integer. So I guess I have two choices: "raster calculator" or "reclassify".

With "raster calculator" I see I obtain my integer raster with the number of classes of thickness I want but the areas are changed a bit (I don't know why).

With "Reclassify" I obtain an integer raster where my areas are essentially the same as before, but it can have a number of classes not higher than the value of thickness it reaches. In other words using the "reclassify" tool  I would have each simulated avalanche with its own scale of thickness and not a common thickness for all my simulations.

Is there a way to accomplish what I'm looking for?


texts: I closed all table but I don't wan't to close underlying map eg like bingHybrid in mapinfo. How to prevent. Any detail about there map and way to handle while closing layers would be appreciated.


texts: Does shapely within function identify inner holes?

texts:
I have a series of shapefiles that include both polygons and multipolygons. I am using ogr to extract them in wkb, and them I am using shapely to check if a point fall inside a polygon (or a multipolygon). However, examining some of the polygons in wkb I realized that many of them have inner rings (i.e. I have a set of coordinates describing the first, outer polygon, and a set of coordinates describing the inner polygon/s). Does the shapely within function take care of this, I mean, does it return True only if a point falls in the area delimited by outer and inner polygons? And are the inner polygons always supposed to be holes?

I am sorry, perhaps these are really dumb questions, but I am completely new to this kind of analyses and, despite my efforts, I could not find some clear information about these doubts. Thanks in advance for your help.


texts: Adding a restriction to network analysis

texts: I currently have some shapefiles regarding flooding that I want to use as restrictions for my network analysis, however, I don't know how to use them. I know I'm able to draw restrictions, but I want to be able to use my shapefiles. Any ideas?


texts: WFS filter for point in polygon (PIP) query

texts:
I am trying to construct a WFS request. I want to determine if a certain point is within a polygon. The point coordinates I want to send via the WFS request. The polygon is saved as a layer on my GeoServer. How does a WFS request like this look like.

My request looks like this, but I am getting an error:

http://myserver.com:8080/geoserver/mapicture/ows?service=WFS&version=2.0.0&request=GetFeature&typeName=mystore:dk_ger_nowater&cql_filter=INTERSECT(geom,%20POINT%20(55.689972%10.239258

The error says:

error on line 4 at column 160: xmlParseCharRef: invalid xmlChar value 16


texts:
I have difficulty to understand how to use "arcpy.SelectLayerByAttribute_management".

I not understand how to write a query.
I want to select the records > 0.

I am attaching my code:

import os,arcpy
shp=(r"C:\Lavoro\\A3.shp")
fil=[f.name for f in arcpy.ListFields(shp)]
print fil
qq=fil[7]
print qq
qq2=fil[8]
print qq2
##exp= '!fil[7]!'>0   #fot test1

##arcpy.SelectLayerByAttribute_management("shp","NEW_SELECTION",fil[7]>0)#test0
##arcpy.SelectLayerByAttribute_management(shp,"NEW_SELECTION",exp)   #test1
##arcpy.SelectLayerByAttribute_management("shp","NEW_SELECTION",'!' + fil[7] + '!'>0)#test2
arcpy.SelectLayerByAttribute_management("shp","NEW_SELECTION",'! + fil[7] + !'>0)#test3

print "End"

I also ask you to explain to me how use   !    "    '  . I don't understand the criteria.


texts: How to delete all rows where a where clause is met with arcpy.da.UpdateCursor

texts:
I would like to use a where clause that finds a certain columns, and if empty then delete the row.

My SQL where clause is;
SELECT UID2 FROM MY311CONTAINERS_2 WHERE SR_SN_2 <> ' ' AND SR_COLOR_2 <> ' ' AND SR_SIZE_2 <> ' ' AND SR_LOCATION2 <> ' '

Here is my code:

import arcpy
aContainerFC2 = 'fc'
fields = ("SR_SN_2", "SR_LOCATION_2", "Dist02", "SR_Color_2", "SR_Size_2")
whereClause = fields + " <>  ' ' "
updCurs = arcpy.UpdateCursor(aContainerFC2, whereClause)
for row in updCurs:
    if not row.getValue(field):
        updCurs.deleteRow(row)

Error:

TypeError: can only concatenate tuple (not "str") to tuple


texts: Layers from a PostgreSQL/ArcSDE database connection are not showing in ArcGIS

texts:
We have a PostgreSQL/PostGIS SDE enabled database, but lately there is a strange issue when adding layers from this database to ArcMap or ArcGIS Pro. We can connect to the database without any issues and see all layers present in it. We can drag any layer in ArcGIS, but the layer is not visible in the data/layout view. It is however visible in the Table of Contents without an exclamation mark. We can even look in the attribute table, no problems there. If I select rows, the contours of the selection are visible, but not the polygon itself. Non registered SDE tables from the same database do work, so I have a feeling there is something wrong with ArcSDE. All tables also work in QGIS.

Is there anything I can do? I could not find similar issues online. Maybe there is something corrupt in the SDE registration?

On the pictures above: left how it is now, right picture is how it should be.

EDIT: I am using ArcGIS Desktop 10.5 and ArcGIS Pro 2.1. Postgres is on version 9.4. It was working just fine before, but suddenly it has this strange bug.


texts: Seeking India population by village, town and city for 1981, 1991, 2001 and 2011?

texts:
Trying to source population data for every village, town and city in India.  I have been able to extract data for 2001 & 2011 from the census site, but am having a tough time finding a reliable source for data from 1981 and 1991.

Any thoughts as to where this may be found?


texts: Rasterio: Use case for src.meta vs. src.profile when reading and writing rasters?

texts:
When I read a raster in with rasterio (version 0.36.0) in Python, I can call up the metadata in two ways:

src = rasterio.open(my_raster)

src.profile

{'affine': Affine(1.0, 0.0, 401900.0, 0.0, -1.0, 7620200.0),
'count': 1,  'crs': CRS({'init': 'epsg:32606'}),  'driver': 'GTiff',
'dtype': 'float32',  'height': 11100,  'interleave': 'band',
'nodata': -9999.0,  'tiled': False,
'transform': (401900.0, 1.0, 0.0, 7620200.0, 0.0, -1.0),  'width': 13750}

src.meta

{'affine': Affine(1.0, 0.0, 401900.0, 0.0, -1.0, 7620200.0),
'count': 1,  'crs': CRS({'init': 'epsg:32606'}),  'driver': 'GTiff',
'dtype': 'float32',  'height': 11100,  'nodata': -9999.0,
'transform': (401900.0, 1.0, 0.0, 7620200.0, 0.0, -1.0),  'width': 13750}

The results are quite similar, but profile stores a little more information such as the tiling and interleave information. Which of these methods is a better practice for retrieving metadata? Are there examples where one is better than the other? Both 'meta' and 'profile' seem to work fine right now for writing the single band float32 raster I am working with.


texts: Error in reproject vector file using GeoPandas

texts:
I have a vector file (in this case shp file) and I read it with GeoPandas (gpd).

vector_folder = r'C:\poly'
vector_file = 'poly.shp'
df = gpd.read_file(os.path.join(vector_folder,vector_file))
df

    id  geometry
0   1   POLYGON ((169628.5983266232 620402.5901412303,...
1   2   POLYGON ((169533.9289514195 620036.5352237758,...
2   3   POLYGON ((169066.8933670809 619834.5738900078,...

The crs is:

df.crs
{'init': 'epsg:2039'}

I would like to reproject it to espg=32636. According to GeoPandas page this should work:

df = df.to_crs({'init': 'epsg:32636'})
df = df.to_crs(epsg=32636)

However, I get this error:

RuntimeError                              Traceback (most recent call last)
<ipython-input-39-1fa06f642778> in <module>
----> 1 df = df.to_crs({'init': 'epsg:32636'})

~\Anaconda2\envs\python3\lib\site-packages\geopandas\geodataframe.py in to_crs(self, crs, epsg, inplace)
    441         else:
    442             df = self.copy()
--> 443         geom = df.geometry.to_crs(crs=crs, epsg=epsg)
    444         df.geometry = geom
    445         df.crs = geom.crs

~\Anaconda2\envs\python3\lib\site-packages\geopandas\geoseries.py in to_crs(self, crs, epsg)
    302             except TypeError:
    303                 raise TypeError('Must set either crs or epsg for output.')
--> 304         proj_in = pyproj.Proj(self.crs, preserve_units=True)
    305         proj_out = pyproj.Proj(crs, preserve_units=True)
    306         project = partial(pyproj.transform, proj_in, proj_out)

~\Anaconda2\envs\python3\lib\site-packages\pyproj\__init__.py in __new__(self, projparams, preserve_units, **kwargs)
    360         # on case-insensitive filesystems).
    361         projstring = projstring.replace('EPSG','epsg')
--> 362         return _proj.Proj.__new__(self, projstring)
    363
    364     def __call__(self, *args, **kw):

_proj.pyx in _proj.Proj.__cinit__()

RuntimeError: b'no arguments in initialization list'

I also tried:

epsg_code = 32636
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg_code)
new = srs.ExportToProj4()
new

'+proj=utm +zone=36 +datum=WGS84 +units=m +no_defs '

df.to_crs(crs=new)

But I get the same error. I get the same error also with GeoJSON formats.

System specifications: anaconda, Python 3, Windows 10, geopandas=0.4.0, proj4=5.2.0


texts: Using gdalwarp to change raster extent, float value issue (Python 2.7, GDAL)

texts:
I am trying to change the extent of a rasterized shapefile to match that of another raster so it can act as a mask for the latter. Other posts suggest using gdalwarp. However, when using gdalwarp to create a new output it throws me a TypeError: argument of type 'float' is not iterable. Is gdalwarp unable to use float values? Although I can't see this written anywhere as a limitation of the function. Below is my code:

slope_src = r'C:\data\extraction\slope.tif'
mask_src = r'C:\data\buildings\blds.tif'
mask_new = r'C:\data\buildings\blds_mask.tif'
slope_raster = gdal.Open(slope_src)
ulx, xres, xskew, uly, yskew, yres = slope_raster.GetGeoTransform()
xres = 0.4
yres = 0.4
lrx = ulx + (slope_raster.RasterXSize * xres)
lry = uly + (slope_raster.RasterYSize * yres)
subprocess.call(["C:\\PRogram Files\\GRASS GIS 7.3\\extrabin\\gdalwarp.exe", "-te", lrx, lry, ulx, uly, mask_src, mask_new])


texts: How to split a linestring every 1000 points in postgis

texts:
To assume compatibility with my survey software, I would like to lighten geometries of my linestrings by cutting them every 1000 points without creating or interpolating points.

I transformed a function found on the net (resample) but it cuts only the first 1000 points...

Does someone knows how to loop it correctly?

CREATE OR REPLACE FUNCTION simplify_npoints(inGeom geometry, maxPoints
integer) RETURNS geometry AS $$

DECLARE
  nPoints integer;
  outGeom geometry;
  fraction float;
  points geometry[];

BEGIN

nPoints := ST_NPoints(inGeom);
  outGeom := inGeom;

IF nPoints > maxPoints THEN
      FOR i IN 1..maxPoints LOOP
      points := array_append(points, ST_PointN(inGeom,i));
    END LOOP;
    outGeom := ST_MakeLine(points);

END IF;

RETURN outGeom;

END;
$$ LANGUAGE plpgsql IMMUTABLE;


texts: Is it possible to create a ternary (triangular) graph template in ArcGIS 10.3?

texts:
Is is possible to edit a .tee graphing template file to create a ternary diagram (example can be found at http://csmres.jmu.edu/geollab/Fichter/SedRx/readternary.html)?

I'm aware of the clever method found at https://www.arcgis.com/home/item.html?id=661a8e7c463a4bd2b529f01221efa8f2. However, that creates a series of shapefiles in a separate data frame. If possible, I'd like to use the built-in graphing tools via ArcGIS.


texts: Improving ArcGIS Results saving and manipulating

texts:
The ArcGIS Results Window is an incredibly useful way to keep track of what has been done in a particular map. I use it constantly as a means of documenting the essential parts of a project's workflow, re-using often accessed geoprocessing tools, and archeology (now, what the heck did we use for the input last time? today's result doesn't $%^#$%& match!).

However, it's also frustrating, because some things which ought to be straightforward just aren't, and it's easy to inadvertently lose stuff you want to keep. Does anyone have any solutions for these friction points? Namely:

ArcCatalog results are not saved at all, if you don't consciously save them they're gone baby.
ArcMap results are only kept if you save an .mxd. If I created an .mxd every time I needed to do a few things my document management overhead would go through the roof. (Imagine needing to create a .doc for every email of moderate import.)
nasty addition to previous: if your ArcMap session crashes, say goodbye to all your results since last save.
Many operations never show up in the Results panel (e.g. Context menu >> Data >> Export). As end users, I suspect there's little to nothing which can be done about this, but it's related to keeping track of steps and is a frequent source of irritation. ;-)
You can drag and drop from Results panel to a model, but not to anywhere else. I want to also copy to a personal toolbar, another map, a folder, a toolbox, a Word document, a ... (Wouldn't it be great if you could just plunk a significant geoprocessing step straight into the metadata editor?)
Some results are Favourites that should are kept indefinitely, e.g. this map project will always need to buffer the selection area of interest by 3km. Currently you can set ArcMap to keep results for a few days, months, or forever. The default of 2 weeks is way too short, and forever will get very very messy and ugly.
Results stay with the map, this is good, but they are also confined to that one .mxd. They're in jail; not so good. It's rare for a project of any significance to have only 1 or 2 .mxd's. Sharing the results is laborious. (This is actually just a restatement of a previous point; light from another angle reveals more facets.)

 * there is a "Shared" item in the results list (see bottom of graphic), but it doesn't seem to do anything. You can't drag a Current or Previous Session into it and it doesn't have a context menu.
((Drag from "Previous Session" working as of 10.3.1. As there is little to know discussion of this elsewhere it may have just been a bug with my local install.))

These result things, the record of steps taken, are important. I want to treat them that way. Not as temporary chaff swept up by the night time cleaning staff on the weekends when I'm not around.

This is a rant, but it's not just a rant. I'm convinced there are means and methods of utilizing what Esri already provides in an improved manner. Let's flush them out of the brush, dust 'em off, and dress up for general consumption.


texts: How to divide polygon into multiple zones given several points (Network Paritition)

texts:
I'm new to QGIS. I have a layer containing all the points of interest. I will like to divide a city into different zones (service area) depending on the location of these points so that each zone will be served by each point. How do I go about doing this in QGIS?
Example: https://www.caliper.com/glossary/what-is-network-partitioning.htm


texts: How to assign mean values from different windows sizes to point feature?

texts:
I have points overlaid on an elevation raster as shown below;

The aim is to get mean value of each point computed from 3x3 window size covering its immediate neighbouring pixels.

It will be highly appreciable if I can append values from multiple raster.


texts: Trying to instantiate geoprocessor with .Net vs 2010

texts:
Retrieving the COM class factory for component with CLSID {5374EC4C-1AA2-4829-A811-DE624ECEC23F} failed due to the following error: 80040111 ClassFactory cannot supply requested class (Exception from HRESULT: 0x80040111 (CLASS_E_CLASSNOTAVAILABLE)).

GeoProcessor gp = new GeoProcessor();
gp.SetEnvironmentValue("workspace", PICKSIMPORT);
IGpEnumList gdbs = gp.ListWorkspaces("", "");
string gdb = gdbs.Next();
while (gdb != "") { Console.WriteLine(gdb.ToString());
 gp.SetEnvironmentValue("workspace", gdb);
IGpEnumList fcs = gp.ListFeatureClasses("", "", "");
string fc = fcs.Next();
while (fc != "") { Console.WriteLine(fc.ToString());
fc = fcs.Next();
gdb = gdbs.Next();


texts: Plugins that works with high resolution screens

texts:
I've written a plugin for QGIS (GeoDataFarm) and I've been using QT creator for all ui files of the plugin, which have worked fine for me. However, I just got an email from someone that tested my plugin on high resolution screen 3200x1800 pixels and apparently it does not work.

This is how the screen look on my computer:

This is how it looks on a high resolution screen:

Any suggestions on how to fix this? Is there any good way of "locking the layout" in QT creator?

Is there any way of testing how it looks on different resolutions, my maximum resolution is 1920x1080?


texts: Sorting integer field in ArcGIS Modelbuilder and using smallest value as variable for tool?

texts:
I have to populate a new field using a formula containing the smallest value of an existing field.

As this is a small part of a larger model I am wondering which is the most efficient method.

My thoughts are as follows:

Sort the field ascending and somehow use the top value
or
Use a formula to search for the smallest integer greater than zero and use this value as a variable?

Which method would be most viable with ModelBuilder from ArcGIS 10.3.1 for Desktop?


texts: Dividing a large raster into blocks for iteration in R

texts:
I have very big rasters for processing (#pixels >1,5*10^10) for this project I am working on. Simple functions such as finding the area of a pixels with a particular value, across the entire raster 1) takes a lot of time and more importantly 2) is unsuccessful because of in_memory problems (yes, I have set rasterOptions(); chunksize() and memorysize())

I think the best ways forward is to divide the raster into blocks and then iterate through these blocks using the doParallel and forEach packages.

I took a stab and am still hitting the same error

Error in rgdal::getRasterData(con, offset = offs, region.dim = reg, band = object@data@band) : long vectors not supported yet: memory.c:3451

This is my code

bs<-blockSize(big_raster)
bs$n #16 blocks
bs$nrow #gives row number at which a chunk can be made

library(doParallel)
library(foreach)

UseCores<-detectCores() #detect number of cores
UseCores #24

cl<-makeCluster(UseCores)
registerDoParallel(cl)

length(cl)

x<-list() #to store areas for different chunks as a list
foreach(ras=1:big_raster[bs$nrow,]) %dopar% {
library(raster)
x[]<-((sum(ras[]==3)*(res(ras) [1])^2)/10000)/10^6 #Area in million ha of all pixels of value 3 in big_raster (linear unit is m)
}
stopCluster(cl)


texts: Index for the spatial distribution/segregation on a raster

texts:
Suppose I have two geotiff raster files for the two time periods t=1,2 where each pixel can get any value in the closed interval [0,1].

Can anyone suggest an index I could calculate to show the change in the spatial distribution of values, e.g. in t=1 pixels with values close to 1 (or at least ≠ 0) where more scattered, while in t=2 such pixels became less homogeneously distributed, i.e. inequality increased? I have stumbled upon the segreg plugin for QGIS, which allows to compute an array of segregation and diversity measures on spatial and non spatial modes, but I am not sure if this can be done with raster data.

Also, how can I calculate a metric which gives the change within the raster extent in the average pixel value?

I am using QGIS 2.18 on a Mac.


texts: Legend widget not loading all layers on inital load esri javascript api 3.20

texts: When using the legend widget (esri JAPI 3.20), not all the layers show up when my map loads. I have to zoom a level first then they all show up. Any fixes?


texts: Non-spatial database analysis in ArcGIS - need alternative solution

texts:
I've recently joined an organisation where we use ArcGIS to do a lot of non-spatial analysis on a spatial dataset. Our main dataset is an environmental polygon survey with lots of attributes (in an ESRI file geodatabase), and I'm frequently asked to conduct analyses such as:
- Calculate the summed area of all polygons which meet certain queries (specified by a Where clause across several fields, for example)

My problem is that I don't know very much at all about traditional database development, but it strikes me that all I'm doing is querying the database (across multiple fields) and then summarizing the results. I've been doing this is in ArcGIS and Python using tools search as Summary Statistics or objects such as searchcursors to scroll through my layers, select records and sum them.

But I imagine non-spatial databases must be a lot faster. But which one? Please tell me what platform I should use to conduct routine summaries on a 500k records sized database. MS Access? Excel and pivot tables? SQL? I don't know where to start, but have a feeling that doing this in ArcGIS (v 10, not 10.1) isn't the best way.

I'm not looking for long answers, just a short shove in the right direction to fill in my lack of knowledge about databases in general.


texts:
I am working with the following ten landsat 8 images:

LC81340412016058LGN00
LC81350412016065LGN00
LC81350422016001LGN00
LC81360412015341LGN00
LC81360422016216LGN00
LC81360432016216LGN00
LC81370412015316LGN00
LC81370422015364LGN00
LC81380412016086LGN00
LC81380422016214LGN00

When I opened these images in QGIS with the supporting shapefiles it covered the entire state of Assam (in India).
But After I tried to merge them using QGIS Raster tool the merged image did not cover the entire state. Only a small part of Assam was shown in the QGIS in the merged image.

Before merging

Before merging in QGIS

after merging in QGIS


texts: Scale on map requires correction coefficient. But why?

texts:
I am computing the distance of a map I am getting from ArcGIS' REST services by using a calculator as provided by this answer. In a nutshell, it uses a GeodeticCalculator for distance computation based on latitude and longitude.

I then use the distance and my display width to compute a scale. My display is currently 362*230 (custom display hence the off measurements). When comparing the distances on screen to the distances as given by Google maps, the measurement always overestimated. However, it did so by a factor of roughly 1.2. I not use this as a correction coefficient and am accurate to about 1m in 1 km.

While this solution works (for now), I would like to know where this offset could come from. Below is a bit of code to show how the computation is done. Note that scale is a field of the class this method resides in.

/** Relevant fields */
private GeodeticCalculator calculator = new GeodeticCalculator(CRS.decode("EPSG:3857"));
private static final double CORRECTION = 1.2;

/** More code left out for reasons of brevity. */

public double computeScale(double latMin, double latMax, double lonMin, double lonMax)
{
    calculator.setStartingPosition(new DirectPosition2D(latMin, lonMin));
    calculator.setDestinationPosition(new DirectPosition2D(latMax, lonMax));
    double distance = calculator.getOrthodromicDistance();
    scale = distance / (Calibration.getTableWidth() / 10) / CORRECTION;
    return distance;
}


texts: Specifying output coordinate system in ArcPy

texts:
I'm updating somebody else's old Python scripts, and trying to make sure I understand their logic and desired outputs while streamlining the code. In this particular script, it specifies an output projection in the processing environment settings.

import arcpy

Output_CS = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',...(etc, etc, etc)"
arcpy.env.outputCoordinateSystem = Output_CS

I think this means that all the subsequent geoprocessing operations (e.g., arcpy.Buffer_analysis) output in WGS84 without needing an additional arcpy.Project_management step.

If that is the case, is there any accuracy risk to running a buffer analysis directly to a geographic coordinate system? In other words, does it do the projecting after the buffer's been appropriately sized?

(The original input data is in a projected coordinate system, one of the NAD83 state plane systems.)


texts: Calculate velocity and errors from distance between GPS points

texts:
I'm using a Kalman filter to track the position of a vehicle and for my measurement, I have a GPX file (WGS84 Format) containing information about the latitude, longitude, elevation and the timestamp of each point given by GPS. Using this data, I computed the distance between GPS points (Using Geodesic distance and Vincenty formula) and, since the timestamp information is known, the time difference between the points can be used to calculate the time delta. Since, we now have the distance and the time delta between the points, we can calculate the velocity (= distance between points/time delta) which could then be also used as a measurement input to the Kalman.

What would be the uncertainty and error of the computed velocity from GPS in comparison to velocity computed directly from vehicle sensors (eg: velocity computed from wheel speed)? (Considering the GPS data itself is noisy, Vincenty is accurate to about 1mm, the geodesic method is accurate to round-off and assuming wheel has no slip, wear and ideal road dynamics)

Code

import gpxpy
import pandas as pd
import numpy as np
from geopy.distance import vincenty, geodesic
import matplotlib.pyplot as plt

"Import GPS Data"
with open('my_run_001.gpx') as fh:
    gpx_file = gpxpy.parse(fh)
    segment = gpx_file.tracks[0].segments[0]
    coords = pd.DataFrame([
    {'lat': p.latitude,
     'lon': p.longitude,
     'ele': p.elevation,
     } for p in segment.points])

"Compute delta between timestamps"
times = pd.Series([p.time for p in segment.points], name='time')
dt = np.diff(times.values) / np.timedelta64(1, 's')

"Find distance between points using Vincenty and Geodesic methods"
vx = []
for i in range(len(coords.lat)-1):
        if(i<=2425):
            vincenty_distance = vincenty([coords.lat[i], coords.lon[i]],[coords.lat[i+1], coords.lon[i+1]]).meters
            vx.append(vincenty_distance)
print(vx)

vy = []
for i in range(len(coords.lat)-1):
        if(i<=2425):
            geodesic_distance = geodesic([coords.lat[i], coords.lon[i]],[coords.lat[i+1], coords.lon[i+1]]).meters
            vy.append(geodesic_distance)
print(vy)

"Compute and plot velocity"
velocity = vx/dt
time = [i for i in range(len(dt))]
plt.plot(velocity,time)
plt.xlabel('time')
plt.ylabel('velocity')
plt.title('Plot of Velocity vs Time')
plt.show()

Reference for GPX Data:
https://github.com/stevenvandorpe/testdata/blob/master/gps_coordinates/gpx/my_run_001.gpx

Plot of Velocity versus Time


texts:
IF there are one or two urban cells in the Moore neighbourhood AND there is a road running through this cell.

How do I add this rule in my raster calculater?


texts: Iterate Rasters will not connect to the Raster Calcuator tool

texts:
Using ArcGIS Desktop, 10.5.1 I need to iterate through hundreds of FGDB rasters (single band, 16-bit, unsigned) that are stored in a single FGDB.  I need to modify each iterated raster with the Raster Calculator and output the resulting raster to a different FGDB.

The problem: I cannot "connect" the Iterate Rasters output to the Raster Calculator tool as an input.  As shown in the image below, there is no input option:

Weirdly, I got the model to work if:

I add the Raster Calculator tool without attempting to connect
I input the Raster Calculator values, as shown below
I hardcode the source FGDB path, as highlighted below (without the hardcoded path, the dreaded 000539 error occurs!?!?!?)

Arc then automatically connects the Iterator to the Raster Calculator, but through the Name bubble.  I thought that the connection always had to be from the "green bubble"!?

Why can't I connect the iterator to the Raster Calculator as an input?


texts: Digitising - Separate fill and line opacities (or remove fill)

texts:
I have a client who specialises in aerial photograph interpretation and is in the process of moving from AutoCAD to QGIS. They find the auto 'fill' when digitising polygons really disruptive and although we can tweak the opacity (see "add feature" tool QGIS 2.0 - polygon fill blocks view of what I want to digitize) ideally we would want a solid line with no fill while digitising.

So far as I am aware there isn't a plugin or other existing solution for this, so my question is how complicated would it be to change the settings to separate opacity of line and fill when digitising? Could this be done via a plugin or only through extending the core QGIS functions?


texts:
I have ortho photos of Romania in Stereo 70 and GPS data in WGS84 of the same country : I want to associate the two series of data.

At the moment, I can't find Stereo 70 in the list of CRS : does it exist ?


texts:
I am working on orthoimagery which I got from drone (UAV). I am wondering if it is possible to create NDVI map from it.

I was following the steps as described in this tutorial:

NDVI in ArcGIS

I got stuck on 1:02 min. In image properties, I was supposed to change the "Band" for Red channel from "Band_1" to "Band_5". I don't have "Band_5" - I have only 4 bands. Does it mean that I am not able to calculate NDVI from my image?


texts:
I need to remove the outside area of a polygon based on a numerical value in a python script.  It would be to determine a core habitat for a species that don't live at the edge of their habitat.  I'm at a total loss on how to approach this.

My ultimate goal is to create a tool in ArcGIS pro, but I need to develop the script first.


texts: Convert ellipsoid height to geoid using SpatiaLite

texts: I would like to convert ellipsoid height to geoid (EGM2008) using SpatiaLite. is there any query that we could achieve this?


texts: I am combining many rasters using a "barcode" that describes attributes for each cell. Combine works on integer values: I combine multiple layers using "combine" tool, then create the bar code by creating a new field (text) and concatenating using the "&" to link fields. The next step is to combine a set of barcodes. But since the bar codes are in a field of type text, "combine" tool will not recognize them. If I copy the barcode values into a field of long integer type, my zeroes create problems. 0000010 becomes 10. I could "trick" the integer field by substituting some other integer for zero. It would perhaps be confusing for the user but it might work. Is there another tool that will combine rasters and keep fields of type text?


texts:
Fairly new to MapInfo both versions. I am a efficient user of ArcGIS, but my query being I was reading up on MapBasic but all the examples for these scripts are for running within MapInfo. Is there any way you can write script (stand alone) outside of MapInfo to make it do something?

I would like some suggestions on what other software can control MapInfo. I am used to writing python scripts outside of ESRI ArcGIS desktop and being able to run the mxds without even opening up the maps. So similarly I wanted to know what can be done for MapInfo?


texts: Elevation profile using linear shape in ArcGIS Desktop?

texts: I'm working with ArcMap 10.2. So, it's not a duplicate of "Displaying elevation profile and storing DEM elevation back to KML line using QGIS?". I have a linear shapefile that has only one element (a river). I have also a DEM. I would like to use my shapefile to make an elevation profile. As it is very long and detailed I really don't want to re-draw it. Suggestions?


texts: Extracting numbers from string field using Python Parser of ArcMap Field Calculator?

texts:
I am trying to calculate field test2 with only the numeric values from a string field test. The entire script is like:

def makestr(test)
     list=[]
     for s in test:
         if s.isdigit():
             list.append(s)
     for a in list:
         str=''.join(list)

But I keep getting an error:

Any ideas?


texts: Tiles map server - How to provide the same tiles on lower zoom?

texts:
I have my own map server defined like below. This is one layer, but have 4 layers for other zooms. Unfortunately for zoom lower then 8 I don't have tile images. Is possible to make zoom out lower than 8 work using the same tiles (as for zoom 8) somehow ? If yes, what can I do, where to start ?

var source = new ol.source.XYZ({
    projection: 'EPSG:2180',
    wrapX: true,
    tileGrid: new ol.tilegrid.TileGrid({
        origin: // ...
        resolutions: resolutions,
        tileSize: [256, 256],
    }),
    tileUrlFunction: function(coord){
        if(coord == null) return null;
        // ...
    }
});

var map = new ol.Map({
    layers: [new ol.layer.Tile({
        source: source
    }), ...],
    target: 'map',
    view: new ol.View({
        center: ol.proj.transform([20.9207687,52.2328546], 'EPSG:4326', 'EPSG:3857'),
        zoom: 8,
        maxZoom: 14,
        // minZoom: 8
    })
});


texts:
After reading replies regarding starting a mapping business, I think I'd better read a Python book.  Any intro-level Python books to recommend?

Mila


texts: Running Standalone PyQGIS Script with Multiprocessing?

texts:
Is it possible to run a standalone PyQGIS script using multiprocessing for better performance?

I'm creating 1,800 rendered map images and would like to run it across 4 processes to speed it up.

If I try to instantiate 4 QgsApplication objects using the following code (one in each process) - they fail silently (i.e. nothing happens and no CPU/hard drive/memory is being used).

QgsApplication.setPrefixPath("/Applications/QGIS3.8.app/Contents/MacOS", True)
app = QgsApplication([], False)
app.initQgis()

If I instantiate it once and then run the 4 processes - I get a "QgsApplication couldn't be pickled" error.

The code uses Python's multiprocessing module to create the parallel processes; using pool.imap_unordered

I'm using QGIS 3.8.2 and Python 3.7.3 on a Mac with MacOS 10.14.6 (Mojave).


texts: What is the unit used by Gdal.GetGeoTransform?

texts:
I tried to retrieve the resolution of a raster image with Gdal in a python notebook to convert surfaces from pixels to hectares.
import gdal
from osgeo import osr

src = gdal.Open(map_raster)
proj = osr.SpatialReference(wkt=src.GetProjection())
_, resx, _, _, _, resy  = src.GetGeoTransform()
src = None

I expected resx to be in meter but a print(resx) gave me 0.00025000000000000017 which would be very small if expressed in meters.
So my question is how can I know the unit system used by the GetGeoTransform function ?
EDIT
the proj value is :
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]


texts: I have downloaded Aura's Ozone Monitoring Instrument (OMI) data from http://mirador.gsfc.nasa.gov website but I can not work with them. How can I open these data in ENVI and what pre-processes are needed for these images?


texts: Loading multiple raster files from list into ArcGIS Desktop?

texts:
How to select and load multiple rasters from one folder in Arc?

My situation is probably familiar to some of you. I do have 3 folders with houndreds of maps in it. I do have also overview sheet, which helps me generate list of rasters that I need. How do I select them automatically from houndreds of others?

Its an issue which occurs for me from time to time. Manual selection is really inconvenient.


texts:
I need to obtain street intersections from reverse geocoded lat/long, where there is no house number information.

I have a tile server and Nominatim running in Ubuntu 12.04 virtualbox.

I followed intructions to install Nominatim v2 at http://wiki.openstreetmap.org/wiki/Nominatim/Installation

Is there a way to modify anything to make this work?


texts: QGIS 3.12 issues unknown CRS on vector shapefile layer

texts:
I'm trying to export in GeoPDF my map from QGIS 3.12.
After updating QGIS 3.6 to 3.12 can't recognize the project CRS EPSG 3120 for raster file and vector layer which is the color layer then I tried to reproject both of layers.
The raster was reprojected correctly with the projection EPSG 3120 but the vector layers still with Unknown CSR :

(Unknown CRS: PROJCRS["Pulkovo_1942_58_Poland_zone_I",BASEGEOGCRS["Pulkovo 1942(58)",DATUM["Pulkovo 1942(58)",ELLIPSOID["Krassowsky 1940",6378245,298.3,LENGTHUNIT["metre",1]],ID["EPSG",6179]],PRIMEM["Greenwich",0,ANGLEUNIT["Degree",0.0174532925199433]]],CONVERSION["unnamed",METHOD["Oblique Stereographic",ID["EPSG",9809]],PARAMETER["Latitude of natural origin",50.625,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",21.0833333333333,ANGLEUNIT["Degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9998,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",4637000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",5467000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER1,LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]] - Projected)

Should I reinstall an older version of QGIS ?


texts: How to get the column header of attribute table

texts: I have a simple vb.net code which gets all the layer names present in the ArcMap and displays it in a vb form combobox. My question is if I select a layer name that is present in my combobox it should get me the column headings of attribute table of the selected layer. How can I achieve this?


texts:
ArcGIS 10.2.2 on Windows 7 32 bit (new Win 7 computer arriving tomorrow).

I'm a complete novice to ArcGIS starting now (June 2014), so speak loud and slow...

I have 10 years of experience with GIS - just not ESRI, so I know what I want to see, now how do I see it??  I've checked the Help, but following the guidance doesn't seem to work.

I'm drawing a small cemetery (4,000 spaces on a few acres).  I want lot labels to appear on the grave lots when I'm at a zoom scale of 1:1 to roughly 1:500 but not any higher.  As the user zooms in and out, these labels should remain one size relative to the polygon in which they pertain.  So when zoomed out they should be rendered in an illegibly small font.  Above that zoom scale (1:501 and further zoomed out) I would like a block number to appear (blocks are on a different layer).  I would like the block number to appear only out to about 1:10,000 and then disappear at 1:10,001 and further out.  At 1:10,001 and further I would like to see section labels appear.

I've been using the standard Label Manager.  Click on Scale Range... and make the appropriate changes to Out beyond and In beyond, but that seems to do nothing.  What am I missing to make this work?

There seems to be a link between these zoom scales and the Frame Reference Scale.  I have tried a Frame Reference Scale of , 1:300, and 1:24,000.  There is no change from one to the other that I can see when I change these.

Also will all these scale zoom visibility features make the jump to a Geodatabase file when I'm finished?  This is the most critical part.  It will be going to an ArcGIS Online account, and this zoom rendering is what I'm shooting for.


texts: QGIS, interpolate along an irregular shape or line to create a DEM

texts:
Interpolating irregular shape, merging DEM's using ArcGIS Desktop asks my question but not for QGIS. Furthermore, the answer is vague and not entirely useful.

My problem is simple. I have xyz data for a river. The river is curved, and I would like to produce a bathymetric map of the bottom. Interpolation is the (nearly) perfect tool. I struggle to get maps that do not highlight the actual data points. The blending is quite poor, and I think this is because the interpolation is happening on a square, or a polygon, depending on which method of interpolation in QGIS you use. The tool does not know that the area of interest is a strange and specific shape!

Is there any way or plugin to tell QGIS to focus the interpolation from a line (which I would draw down the center of the river, or a polygon, (which I could draw outlining the river) in order to interpolate in a very specific area?

I do not want to use a polygon to clip my interpolation result. It is extremely close to what I want but I feel like the rest of the interpolation that I do not want is affecting the river.

I have a feeling GRASS can do this. If someone can confirm this and tell me the function I will learn how to use GRASS (although I'm kind of scared of it). Perhaps there is a function or plugin in QGIS.

Adding this image to help with a comment I made. Demonstrates the problem and also another issue with the actual path of the sensor being too visible.


texts: Creating ArcGIS tool that can run list of checks?

texts:
Is it possible to create a tool in ArcMap where by uploading a shapefile, the tool can run a various number of checks/tests (ie intersecting other shapefiles) and flags or highlights conflicting datasets upon initiating?

Maybe even giving a further option of clipping of the shapefile around the flagged dataset?

I’m a student and I’ve found that this might be possible with arcpy, however I don’t exactly know the limitations of using this method and even if such a tool could be made.


texts: Splitting a line shapefile into segments of equal length in python

texts: I want to split a line shapefile into segments of equal length, let say into length of x meters each. There is tool v.split.length as mentioned in this answer: how to split a line into a specific no. of parts? Is there a way to do that in python without using QGIS? May be using shapely and fiona? as there is an interpolate method in shapely which could be used to add points at a specified distance on a line. But I couldn't find a way to divide a line into segments of specified length.


texts: Creating a new field based on ranges in a field using Field Calculator in QGIS

texts:
I have a field (Tot_Score) that has scores from 2,3,4,5,6,7,8. I would like to create 4 classes out of these scores whereby the classes will be created in another field (Vul_Class) so 1st class = 2 and 3, 2nd Class = 4 and 5, 3rd Class = 6 and 7, 4th Class = 8.
What is the script I would need in the QGIS Field Calculator to do this?


texts: Auto resize map to 100% after content pane is removed in ArcGIS JSAPI

texts:
Been at this for hours and I can't seem to figure it out.
I'm using AG-JSAPI 3.6

I have a map layout with a bordercontainer and two contentpanes, One is the map div and the other is the left div (which holds the legend, and basemaps and such as child elements).

Anyway, I have a checkbox element on the footer that when clicked, fires off a function that (should ideally) make the left div disappear and the map div expand to the width of the browser.

function toggleDiv() {
    // If statement to check if the checkbox is checked.
if (dijit.byId("ToggleLeft").checked) {
    //Sets the left div's display CSS property to inline and then setting the
    //width property
   dojo.setStyle(dijit.byId("leftDiv").domNode, 'display', 'inline');
   dojo.byId("leftDiv").style.width = '285px';

}
else
{
    // Sets the display css property to "none" which makes the element disappear
   dojo.setStyle(dijit.byId("leftDiv").domNode, 'display', 'none');

};

}

The thing is, when I resize the browser, the map automatically resizes to the width of the browser. However, I want to do this as soon as the leftDiv element disappears. What js code would I have to use to get this to happen? Any help would be appreciated.

Thanks


texts:
I'm trying to project from ED50 31N to WGS 84. I'm using script below but when I generate the outputs they are inverted from the coordinates I would expect. I expect x_out to be 4.031819616199967 and y_out to be 58.63423101793458 and not vice-versa (what python returns for me). I'm running this with pyproj 2.6.1. When I use an online converter (e.g. https://epsg.io/transform#s_srs=23031&t_srs=4326&x=560000.0000000&y=6500000.0000000) I do get the expected coordinates. Any idea why they are inverted with this script?
Python 3.8.2 (default, Apr 27 2020, 15:53:34)
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license()" for more information.
>>> from pyproj import Transformer
>>> in_EPSG = 'epsg:23031'
>>> out_EPSG = 'epsg:4326'
>>> transformer = Transformer.from_crs(in_EPSG, out_EPSG)

>>> x_in = 560000
>>> y_in = 6500000
>>> x_out,y_out = transformer.transform(x_in,y_in)

>>> print(x_out)
58.63423101793458
>>> print(y_out)
4.031819616199967


texts: OpenLayers crashes stock browser on Android 4.0.4 devices

texts:
Update: The issue appears to be with OpenLayers Release 2.13 dev. I mistakenly mentioned 2.12 as being problematic but have retested with no problems. Also, no problems with 2.13rc3

Update 2: This issue appears to exist again in 2.13.1

We're developing a mapping tool using Openlayers 2.1.2 (issue also confirmed on 2.1.3 ).

I've tested on 4 devices running Android 4.0.4. The map panning/zooming works for a few seconds before the browser completely locks up - no longer able to pan/zoom or navigate to another page. The browser then either crashes or I need to force stop in application settings.

Devices running 4.1.1 have no problems.

This is also an issue with OpenLayers own mobile example page:
http://openlayers.org/dev/examples/mobile.html

Again, after a few seconds, the browser freezes.

I see no reports of this online. The OpenLayers documentation says it is supported on Android 2.1+

Has anyone else run into this issue?


texts: Changing layer symbol across multiple scales in ArcMap?

texts:
I have some Cycling paths as line features, and I need to adjust the line visibility. When zoomed 1-1000 I need to see all lines. When zoomed 1000+ I need to see only red lines. Of course not only colors but attributes are various. So in that adjustments I can view exact path contours, but zoomed out 1000+ I only see the location.

Red lines are center lines of the blue lines.

I can't export to other layer, because I need to upload layer to a server, where is one layer accepted.


texts: Create prj-file for shapefile in QGIS 3?

texts:
When you load in QGIS a shapefile without a PRJ file, you can set the coordinate system, but no PRJ file is created.

This older question: How to obtain .prj files in QGIS?

has an answer saying menu "Vector - Data Management Tools - Define Current Projection" does it, but it is not present in recent versions. Any suggestion of how to do it?

Of course, I can do it manually, copying the PRJ file of some other shape, but probably there is a way of doing it from within QGIS menus.


texts: Removing a point feature that seems protected against removal

texts:
What if changes to a certain dataset are protected?

Can't remove a faulty point in the regular way. The dataset is created by myself in the field. A set of GPS locations imported as a KML file. One is just offline the transect.


texts: Class Probability Tool in ArcGIS Desktop 10.2

texts: I am trying to create a multi-band raster using the class probability tool in ArGIS Desktop 10.2. I used the classified image I created based on my signature file as the input raster, and my training classes (.gsg) extention as my input signature file. All I get as an output when the tool runs is a mostly grey layer with no actual information attached to it. Anyone have any ideas about what I am doing wrong.


texts: How does Nominatim know my address which is not (as far as I can tell) in OSM Database?

texts:
I am looking for an additional source for geocoding and reverse geocoding.  I went to http://www.openstreetmap.org and clicked the button to query my current location.  Nominatim returned a "house" with my correct address.  When I search for my address, I get the same point.  So I figured that someone must have entered a node with building attribute of house and my address.  But when I download the OSM data for my state into QGIS, and then create a point layer with it, when I filter the layer on nodes with "building"="house" set, there are only 330 or so in the entire state, all spread out.  Like only people who have added their own house have these nodes, perhaps.  The majority of nodes that have "addr:streetnumber" set are for businesses, churches, and other POI entries.

Which brings me back to the original question - how does Nominatim know of my house?  FYI I also created a layer of the closed polygons, thinking the house may be there, but it is not there, either.


texts: Replacing ID representing field with actual values depicted by ID when styling with data defined properties or Labeling?

texts:
I am currently using data defined properties and labeling to present pipes based on diameter.
I now have a dictionary table that has id's for every combination of diameter and material and I'm using it with the value relation widget.

I must find a way to access the value represented by the key in the diameter || material field.

An ugly way that solves this would be to use:

CASE
   WHEN id_material_diam in (id1, id2) THEN result_formula_containing_diameter
   WHEN id_material_diam diameter in (id3...) THEN ...
END

Is there an elegant way that I can access the actual values for which the ID is the key?

I am using QGIS 2.2 Valmiera


texts: Making buffer zone based on size of attribute using ArcGIS for Desktop?

texts:
Buffer zone isn't really the correct word, sorry, but I am not sure how to explain it.

I want to make a circle around each attribute based of the size of the value, similar to a bubble map. However, I also need to select the overlaps to see which area has the largest overlap (does this make sense?)

I am using Vector data (I think) in the form of xlsx datafile (need to geo-code this first) and I do have access to a large amount of GIS programs (school computer)

How should I go about doing this?


texts: Changing data source of layers with ArcPy?

texts: How do I change the data source of all layers in a map document using ArcPy and Python?


texts:
I'm trying to save a custom layer by doing the following

ILayerFile layerFile = new LayerFileClass();
layerFile.New(saveFileDialog.FileName);
layerFile.ReplaceContents((ILayer)customProperty);
layerFile.save();

The custom property is actually CustomLayer2 which is nothing but a class that extends the BaseCustomLayer

It gets to layerFile.save() but then craps out.  I get the following error

  Error: HRESULT E_FAIL has been
  returned from a call to a COM
  component

I'm not sure if this is because I'm trying to save a custom layer or not.  But since BaseCustomLayer extends the ILayer shouldn't this be all good?  It's not failing at ReplaceContents and it looks like data is actually getting there.  The only field that's null is ThumbNail.  Can this be null?  Does it need to be created in the CustomLayer class?


texts: QGIS OpenFileGDB compatibility with File Geodatabase?

texts:
I am trying to view an ESRI FileGDB in QGIS 2.8.1.

Do I have to use "add vector layer" - "directory" - "OpenFileGDB" in order to open and view the data?

Is there a way to specifically open the ESRI FileGDB type in QGIS 2.8.1?

The data I am trying to view is not being displayed correctly.  It appears different in other applications.  The polylines appear generalized even though I turned off simplify geometry.
Two layers in which the polylines should be identical are being displayed differently.  The polylines do not match and appear to have different vertices making up the polylines.

I got shapefiles of the same data from the same source and they appear correct.
The projection is set to the correct projection in QGIS.

Any ideas on how I can get the ESRI FileGDB to display correctly in QGIS 2.8.1?


texts: geopandas' .crs doesn't store crs on 'proj' file

texts:
As occurred in this question: Add CRS to a shapefile in Geopandas?

f.shp has not crs

This code:

 shp=gpd.read_file('f.shp')
 shp.crs = {'init' :'epsg:2100'}
 out=path
 shp.to_file(out)

doesn't store the crs in the proj file and when I open it in Arcmap it says that it doesn't have a crs.

Although

shp.crs

gives

{'init': 'epsg:2100'}

How can I make it have the assigned crs when I export it?

I also tried:

shp.to_file(out,crs_wkt='2100')

so it provides the crs while exporting but that didn't fix it.


texts: while running fill sinks in QGIS 2.14.6 i am getting this error..what is the reason?

texts:
Algorithm Fill sinks starting...
io_gdal 0 -TRANSFORM -INTERPOL 0 -GRIDS "C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\c8603b46118c467c8c9d349dd289bde8\dem11.sgrd" -FILES "F:\my projct\dem11.tif"
ta_preprocessor "Fill Sinks (Planchon/Darboux, 2001)" -DEM "C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\c8603b46118c467c8c9d349dd289bde8\dem11.sgrd" -MINSLOPE 1 -RESULT "C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\0142c11dfb5940c2b88c132dfe44efdd\RESULT.sdat"

C:\PROGRA~1\QGIS2~1.14\bin>set SAGA=C:/PROGRA~1/QGIS2~1.14/apps\saga

C:\PROGRA~1\QGIS2~1.14\bin>set SAGA_MLB=C:/PROGRA~1/QGIS2~1.14/apps\saga\modules

C:\PROGRA~1\QGIS2~1.14\bin>PATH=c:\progra~1\qgis2~1.14\apps\python27\lib\site-packages\shapely-1.2.18-py2.7-win32.egg\shapely\DLLs;C:\PROGRA~1\QGIS2~1.14\apps\Python27\DLLs;C:\PROGRA~1\QGIS2~1.14\apps\qgis-ltr\bin;C:\PROGRA~1\QGIS2~1.14\bin;C:\Windows\system32;C:\Windows;C:\Windows\WBem;C:\PROGRA~1\QGIS2~1.14\apps\Python27\Scripts;C:/PROGRA~1/QGIS2~1.14/apps\saga;C:/PROGRA~1/QGIS2~1.14/apps\saga\modules

C:\PROGRA~1\QGIS2~1.14\bin>saga_cmd io_gdal 0 -TRANSFORM -INTERPOL 0 -GRIDS "C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\c8603b46118c467c8c9d349dd289bde8\dem11.sgrd" -FILES "F:\my projct\dem11.tif"
ERROR 1: Can't load requested DLL: C:\PROGRA~1\QGIS2~1.14\apps\saga\dll\gdal_MrSID.dll
127: The specified procedure could not be found.

ERROR 1: Can't load requested DLL: C:\PROGRA~1\QGIS2~1.14\apps\saga\dll\gdal_MrSID.dll
127: The specified procedure could not be found.

## #####

### ##

# ## ## #### #

##### ## #

# ## ##### #

library path: C:\PROGRA~1\QGIS2~1.14\apps\saga\modules\io_gdal.dll
library name: Import/Export - GDAL/OGR
tool name : GDAL: Import Raster
author : O.Conrad (c) 2007 (A.Ringeler)

Parameters

Grids: No objects
Files: "F:\my projct\dem11.tif"
Transformation: yes
Interpolation: Nearest Neighbor

loading: F:\my projct\dem11.tif

Driver: GTiff

Bands: 1

Rows: 4789

Columns: 4196

Transformation:

x' = 393016.204524 + x * 30.555556 + y * 0.000000

y' = 2596593.474508 + x * 0.000000 + y * -30.555556

warning: top-to-bottom and left-to-right cell sizes differ.
using cellsize: 30.555556

loading band [1/1]

C:\PROGRA~1\QGIS2~1.14\bin>saga_cmd ta_preprocessor "Fill Sinks (Planchon/Darboux, 2001)" -DEM "C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\c8603b46118c467c8c9d349dd289bde8\dem11.sgrd" -MINSLOPE 1 -RESULT "C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\0142c11dfb5940c2b88c132dfe44efdd\RESULT.sdat"
ERROR 1: Can't load requested DLL: C:\PROGRA~1\QGIS2~1.14\apps\saga\dll\gdal_MrSID.dll
127: The specified procedure could not be found.

ERROR 1: Can't load requested DLL: C:\PROGRA~1\QGIS2~1.14\apps\saga\dll\gdal_MrSID.dll
127: The specified procedure could not be found.

## #####

### ##

# ## ## #### #

##### ## #

# ## ##### #

library path: C:\PROGRA~1\QGIS2~1.14\apps\saga\modules\ta_preprocessor.dll
library name: Terrain Analysis - Preprocessing
tool name : Fill Sinks (Planchon/Darboux, 2001)
author : Copyrights (c) 2003 by Volker Wichmann

Load grid: C:\Users\WATERR~1\AppData\Local\Temp\processing1d1d4cddaf354226b0924459ec8556cf\c8603b46118c467c8c9d349dd289bde8\dem11.sgrd...

Parameters

Grid system: 30.555556; 4790x 4197y; 393031.482302x 2468367.085619y
DEM: dem11
Filled DEM: Filled DEM
Minimum Slope [Degree]: 1.000000

C:\PROGRA~1\QGIS2~1.14\bin>exit
Converting outputs
Loading resulting layers

The following layers were not correctly generated.
Filled DEM
You can check the log messages to find more information about the execution of the algorithm


texts: Geoserver Raster layer cropped by polygon?

texts:
I've got a Geoserver raster layer which represents a satellite image. I use the Geoserver WMS to display the raster in OpenLayers. I want to display only parts of the raster - those within a certain region.

Is it possible to send a parameter to the server (like a bbox, or a CQL expression), which will tell it to render only within a polygon?

Can i make any 'tweaks' in OpenLayers to accomplish this behavior?
For vector layers I can use a cross-layer CQL filter in the WMS request, but what to do with a raster?


texts: Understanding odd behaviour of QGIS symbol picker?

texts:
I am working on a symbol set for Swedish road signs. I have batch converted some 300 png files to svg, and they work very well in QGIS.

Except for one detail: Only the top half of the symbol is visible in the symbol picker (the signs below are actually triangular and round). Looking at the code, these files are no different from other svg files in the symbol library. It makes no difference if I scale down the images.

Is there some svg parameter I need to add?

I used Imagemagick for the batch conversion, and QGIS is 2.18.4.


texts: Copy parts from several attribute tables in shapefiles?

texts:
I have several shapefiles from which I wish to take parts from their attributes tables to make a new attributes table.

Is this possible without destroying the original attributes table? I have features connected to them.

If I were to use the join tool, how do I join from several layers?

How can I do this?


texts: Connection issue with ArcGIS Server 10 map services with 931 SDE

texts:
I am trying to figure out an issue with map services. I have non pooled services running on ArcGIS server 10 environment. These services connect to 931 SDE. I am observing that map services (even not in use) are consuming SDE connections and ultimately connections are maxed out. I am looking for some recommendations on using 931 SDE with 10 map services.

Thanks
jay


texts: gdal_translate to JPEG has no EXIF in the output

texts:
I am trying to create .jpg files that can be consumed on the web from a .tif image. I would like to rely on EXIF data about the image metadata, however, I fail to read any. What am I doing wrong?
$ gdal_translate --version
GDAL 3.0.4, released 2020/01/28

$ gdalwarp -t_srs EPSG:4326 IW32001452_mask.tif /tmp/tif2jpeg.tif
Creating output file that is 1802P x 1244L.
Processing IW32001452_mask.tif [1/1] : 0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdal_translate -of JPEG -scale -co write_exif_metadata=yes /tmp/tif2jpeg.tif IW32001452_mask.jpg
Input file size is 1802, 1244
0...10...20...30...40...50...60...70...80...90...100 - done.

$ exif IW32001452_mask.jpg
Corrupt data
The data provided does not follow the specification.
ExifLoader: The data supplied does not seem to contain EXIF data.


texts: Export selected columns (to excel or dbf) in ArcGIS Desktop 10.1

texts:
I have lots of column in shapefile and I need few columns.

Therefore, I am trying to select using control and trying to export but it exports all the columns (I did not get “selected records “option while exporting)

Image 1 :

I tried using ArcToolbox (Data management tools and few more) but not getting suitable option.

Researching on Google and going through following posts:

Post 1 ,Post 2 (I do not think its proper way)

So any help will be great!!


texts:
I am working on two shapefiles that I need to modify in ArcGIS.

Sometimes ArcGIS simply doesn't open the shapefiles. The error shows the opening feature class message.

Then, while trying to save them, if I just use the Save as a command, I get a .xmd file.

Now, I go to the layer and try Export Data, and re-write the output feature class from "\\Mac\Home\Documents\ArcGIS\Default.gdb\..." to "\\Mac\Home\Documents\ArcGIS\blablabla.shp".

I got the shapefiles, however; there was a random one with some polygons missing, etc. and which it doesn't take into account my modifications.

Does someone know how can I fix this?


texts: Buffering in meters/km using WGS84 layers using QGIS

texts:
I am using QGIS 1.7.4-Wroclaw and working with an SHP layer in WGS84. I want to buffer out to a certain number of kilometers from this layer using the ftools buffer tool.

I understand that the buffer tool always uses the layer units, which for WGS84 are decimal degrees. We all know that degrees don't convert consistently to meters, so how should I go about making my buffers?

Is it necessary for me to convert the shapefiles to a different CRS that natively uses meters?

If so, how do I choose one?


texts: How to configure "Same as Display" as processing extent by default in ModelBuilder?

texts:
I'm working on a very simple model for creating a path from a selected point to the closest polyline, in this model I used a Cost Distance and Cost Path tool for creating the more feasible track (considering the slope).

Everything works fine, but because I'm using a Cost Raster from a very large area (whole México country) the processing is very slow, so I set the processing extent as "Same as Display", the problem comes when I select another point from a different extent, the processing extent have to be configured in the model as "Same as Display" other way it is automatically set as "As Specified Below"

What I want to know is if there is a way to configure the Processing Extent as "Same As Display" and it stay like that as default doesn't matter where I select the point. I want this model to be published on web but I don't want to add the processing extent as parameter


texts: I have a polygon feature class representing discontinuous fields (see below), I would like to make a heatmap based on the distribution of these polygons. I'm using Arcmap 10.2.2 and am aware of the methods to make a heatmap from point data but I'd like the heatmap to take into consideration the size of each site.


texts: Using ogr2ogr to export an Oracle Spatial table with measured geometries (3302) to a shapefile ignores the measure values

texts:
I have an Oracle spatial table with measured geometries (SDO_GTYPE=3302) which I want to convert to a shapefile.
The command I am using is this:
ogr2ogr -f "ESRI Shapefile" "tshapeARCM.shp" -lco SHPT=ARCM -dim XYM -nlt MULTILINESTRINGM OCI:{username}/{password}@{db_alias}:{table_name} -progress -skipfailures

I tried to use any combination of the lco/nlt/dim parameters but the shapefile insists on coming out with all M values as 0
I also checked the layer using ogrinfo and it reports:
1: {table_name} (3D Unknown (any))

Does this means that it does not recognize the 3302 geometry? Is this a bug in ogr2ogr or I am missing something?
Using GDAL 2.2.3 on a windows 10 laptop


texts: RuntimeError: cannot open "feature class" with SearchCursor?

texts:
I have a very simple script that isn't working.

import arcpy

workspace = "C:\\Users\\Work\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\master.sde"
arcpy.env.workspace = workspace

wire = "gis.DBO.wire"

with arcpy.da.SearchCursor(wire, "id") as cursor:
    for row in cursor:
        print row[0]

It returns RuntimeError: cannot open 'gis.DBO.wire'

UPDATED

I still haven't figured this out. I just did a repair on ArcGIS Desktop 10.3, and that didn't solve it. So, then I removed Python, and reinstalled it and it still doesn't work.

I have this script I've tried running in PyScripter and IDLE:

import arcpy
sde = "C:\\Temp\\Data.sde"
domains = arcpy.da.ListDomains(sde)

I get:

And, yes it is a valid sde connection.

If I do the exact same code above in the Python window of ArcCatalog it works fine.

Here's the IDLE error:

Traceback (most recent call last):
  File "<pyshell#2>", line 1, in <module>
    domains = arcpy.da.ListDomains(sde)
RuntimeError


texts: draw extend of print composer map as rectangle in canvas

texts:
After defining the extend of my printed map in the printcomposer I can show the map center in the canvas by hitting "view extend in map canvas". So far so good. But is it also possible to show the exact extend with all borders? Or even better: map composer extend is shown as rectangle in the canvas.

At the moment I get that result by copying the coordinates of the map-composer extend to a txt-file (;-sparated), import it to the canvas and connect the points as polygon. But if I change the extend of the print composers map, I have to do all steps again. So is there a other way?


texts: How to convert data in NAD83 (NSRS 2007) to NAD83 (CORS 96)?

texts:
I'm quite a ways outside my box on this one, but I would ceratinly appreciate any perspectives on this issue.  I have some control points from NGS that are in NAD83 (NSRS 2007) that I would like to merge with the rest of our control point dataset, but our existing data are in NAD83(CORS 96) ... I think ... or it could be NAD83 (2001).  Perhaps I should ask:

Which NAD83 is inherent in ESRI's NAD_1983_StatePlane_Georgia_West_FIPS_1002_Feet coordinate system (Geographic: GCS_North_American_1983  Datum:   D_North_American_1983)?

In any case, does anyone know of any transformations for converting NAD83 (NSRS 2007) to one of the "regular" NAD83 datums?  NGS will not be publishing any, and ESRI does not have a transformation filter for it.  Any help would be appreciated.  BTW, I did read the GIS-SE posts from last year regarding NAD83<-->HARN conversions, but I am wondering if anyone has dealt successfully with this issue since then.


texts: How can I fix broken links in an mxd with only an ArcGIS Server license?

texts:
I've been asked to fix a broken arcgis server site.  I suspect the problem is that the paths in the mxd don't point correctly to the data files.

Without installing desktop is there some way to confirm that the paths are correct, and fix them if they are broken?


texts: How to force sql server to use the spatial index?

texts:
I'm beginner with SQL Server (2012). I created a spatial index, and I'm trying to force SQL Server to use that index to solve the query.

Here's what I have:

DECLARE @g geography

SET @g = geography::STPointFromText('POINT(-79.25161 43.735965)', 4326)

SELECT *

FROM CUSTOMER

WITH(INDEX(idx_Geography))

where HomeGeog.STDistance(@g)

The spatial index (idex_Geography) is created for HomeGeog.

I'm getting this error:
Msg 102, Level 15, State 1, Line 5
Incorrect syntax near '('.

I'm not sure where the error is...

Thanks for you help!


texts: Moving polygons to points and rotated using QGIS

texts:
I'm using QGIS 2.18.11. So I have a list of various coordinates and rotation angle. I'd like to paste a template polygon with each of their center point positioned on each coordinate and rotated by an angle as specified in the list.
For the translation part, I tried to use the method outlined here: Move a layer in QGIS to a precise position based on centroid but got stuck near the end because QgsFeature.setGeometry() is voided.

feat = QgsFeature()
feat.setGeomertry(QgsGeometry.fromPolygon(pol_t))
Traceback (most recent call last):
  File "<input>", line 1, in <module>
AttributeError: 'QgsFeature' object has no attribute 'setGeomertry'

As for the rotation part, I came across QgsVector.rotateBy() that seems to be able to do the trick. Am I on the right track?


texts: QGIS Georeferencer - Coordinates Rounding Too Small

texts:
I have an image I'm trying to georeference on my map, but am having an issue with the coordinates which are sent to the map from the georeferencer.

As you can see below in the image, the yellow boxes above each point each show the correct location where the points should be referenced to. An example being 37.2317089 rather than 37.23 as shown in the table in the picture below.

Is there a way to increase the amount of digits after the decimal point?

Or do you have any suggestions for any other ways of going about this?


texts: Simulating Random Agent Movements along Roads on Map?

texts:
I need to simulate agent movements on map as follows:

Agent starts a known location and moves along a road. At each road intersection, the agent uniformly chooses one possible road and moves along that. The path chosen is to be traced with a poly-line.

I know how to choose a path uniformly and mark it with poly-lines once I have the data. But I am looking for how I can retrieve data like road intersection points and the roads leading from there on.


texts: Verifying if user input as value exists in column using ArcPy?

texts:
I am trying to verify whether a value inserted by a customer exists in the field.

import arcpy
import os
import sys

#Set the work space
arcpy.env.workspace = "Y:\Milad\GeometricNetworkAnalysis"
path= arcpy.env.workspace

#http://pro.arcgis.com/en/pro-app/arcpy/get-started/data-access-using-cursors.htm
user_boundary = raw_input("Insert the number of the Node? ")
Table = "RFNodeBoundary.shp"
field = ['Name']

with arcpy.da.SearchCursor(Table, field) as cursor:
    for row in cursor:
        #print (u'{0}'.format(row[0]))   #this is working
        if user_boundary in cursor:
            print "your node is here."
        else:
            print "The node is not in the database."

The code returns "The node is not in the database." whether the node exists in database or not.

How I can fix the problem?


texts:
I have two SQL table with same columns/field name, but when I open both tables in ArcMap it shows different information in Source Tab.

SDE table type shows all of the fields where query table type missing one field in ArcMap though all same field I can see from SQL Server end.

I am using ArcMap 10.2 also I checked in pro and same thing is happening.

I'd like to know

How I can change Query Table to SDE table.
Any idea why two identical tables with field name, one showing all field and
other one is missing which is Date type column?


texts: Best way to correct UAS imagery sun glint over a water surface

texts:
I have a challenge of where there is a sun glint which includes both specular and non-specular  brightening of a water surface on Unmanned Aerial Vehicle (UAS) imagery. Please see the raw image below before vignette correction is applied. The sun glint is the lower orange and red part of the image. The red spots are specular reflection and the orange pixels are brightening due to the sunlight.

I am able to identify the sun glint with the following code, which is not perfect as it requires seeing the image and estimating the percent of the sun glint on the map which is variable. I am dealing with thousands of images, which makes it difficult to use a single percent for all images as level of sun glint varies.

glint_min_value = np.percentile(band.ReadAsArray(), 85.0)
final_gl_tif = glint_tif.replace('_gl', '_fi_gl')

if not os.path.isfile(final_gl_tif):
    cmd = 'gdal_calc.py -A "{}" --quiet --outfile="{}" --calc="A*(A{}{})" --NoDataValue=0'.format(
        img, glint_tif, operator, glint_min_value, nodata
    )
    call(cmd, shell=True)
    print('Glint tiff {} {}'.format(glint_tif, os.path.isfile(glint_tif)))
    # project and set 0 to -9999 nodata value

    gdal_wrap = 'gdalwarp -q -srcnodata {} -dstnodata {} {} {}'.format(
        SOURCE_NODATA, nodata, glint_tif, final_gl_tif)

    call(gdal_wrap, shell=True)

I would like to find away to adjust the image so that the sun glint is effect is corrected. Is there any known method or approach for this?


texts: Concatenating metadata from Features in a FeatureCollection

texts: I have a FeatureCollection made up of many features with single piece of data in the properties like {date:ndvi}. Many of these features have the same geometry(polygon). Is there an easy way to concatenate features with the same geometry into single features that has a single properties dictionary but all metadata? i.e. 5 features with same geometry and properties {date1:ndvi1}, {date2:ndvi2}, {date3:ndvi3}, {date4:ndvi5}, {date5:ndvi5} becomes single feature with properties {date1:ndvi1, date2:ndvi2, ... , date5:ndvi5}.


texts: Data Migration: PostgreSQL (PostGIS) to SQL Server 12

texts:
I'm migrating data from a PostgreSQL database (with PostGIS), which contains geospatial data to SQL Server '12 and having some import errors. I'm assuming these errors are due to the geospatial data types.

I originally ran pg_dump dbname > outfile and then loaded the output file which is  a .bak into SQL Server '12 using the GUI import wizard.

Import Errors:

Error 0xc02020a1
Error 0xc020902a
Error 0xc0202092
Error 0xc0047038

Question: What is the best way to migrate PostgreSQL geospatial data into SQL Server '12?


texts:
I would like to cluster a point cloud and then export the centre and label of each cluster to a text file.
So far I have this as my pipeline:
{
    "pipeline":[
        {
            "type":"readers.rxp"
        },
        {
            "type":"filters.range",
            "limits":"reflectance[0:25]"
        },
        {
            "type":"filters.dbscan",
            "min_points":10,
            "eps":0.05,
            "dimensions":"X,Y,Z"
        },
        {
            "type":"filters.range",
            "limits":"ClusterID[0:]"
        },
        {
            "type":"filters.groupby",
            "dimension":"ClusterID"
        },
        {
            "type":"filters.stats",
            "dimensions":"X,Y,Z"
        },
        {
            "type":"writers.text",
            "format":"csv",
            "order":"X,Y,Z,ClusterID",
            "keep_unspecified":"false",
            "filename":"outputfile.txt"
        }
    ]
}

But this produces a .csv with all cluster points listed.


texts: Export image from Google Earth Engine to my drive

texts:
I  have created code that takes the reflectance of Landsat8 and clip it with Iraq's borders. the problem is that I can't export it to drive and I get this error messege:

  Scale must be greater than 0.010000 m/pixel

I took the scale from the data catalog, but it may be that I haven't fully understood the scale idea/logic.

This is the code I wrote:

// Load a Iraq boundary from a Fusion Table.
var Iraq = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw')
  .filter(ee.Filter.eq('Country', 'Iraq'));
print(Iraq)
Map.addLayer(Iraq, {}, 'Iraq from fusion table');

// Load a Landsat 8 -atmospheric corrected
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
//surface reflectance
    .filterBounds(Iraq)
    .filterDate('2015-03-01', '2018-03-31');

// Get the number of images.
var count = collection.size();
print('Count: ', count);

// Sort by a cloud cover property, get the least cloudy image.
//var image = ee.Image(collection.sort('CLOUD_COVER').first());
//print('Least cloudy image: ', image);

// make a composite (mean image) of the images
var meanImage = collection.mean().clip(Iraq);
print(meanImage);

//make image for mesopotamia only
var mesopotamia=meanImage.clip(table);

//visualization
var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

Map.addLayer(meanImage,visParams,'Landset8');
Map.addLayer(mesopotamia,visParams,'Mesopotamia');

  // Export the image, specifying scale and region.
Export.image.toDrive({
  image: meanImage,
  description: 'IraqLandsat8',
  scale:0.0001,
  region: Iraq
});

  // Export the image, specifying scale and region.
Export.image.toDrive({
  image: mesopotamia,
  description: 'IraqLandsat8',
  scale:0.0001,
  region: table
});


texts: Cause of ArcGIS Raster Calculator Parsing error : invalid syntax?

texts:
This is my first attempt at much code in the raster calculator of ArcMap 10.0, but I think I'm close.

I have nine rasters (A-I below) and want a new raster with value = 1 where any of the rasters have field 1 >= 1 and <= 366 and field 2 not equal to 5.

Here's what I've got:

Con(([A].1|[B].1|[C].1|[D].1|[E].1|[F].1|[G].1|[H].1|[I].1 >= 1  &  <= 366) & ([A].2|[B].2|[C].2|[D].2|[E].2|[F].2|[G].2|[H].2|[I].2 != 5) , 1,0)

I get the little red circle with an X in the top left of the raster calculator window and an error message:

  Parsing error : invalid syntax (line
  1).

Anyone know what I've got wrong?


texts: Calculate and Restart Sequential ID for Unique Values in a Field

texts:
I need to calculate a sequential ID for unique values in a field. So, for the values in Field1, the sequential ID count should count all occurrences of that value, then restart the count for the next value. For the example below, column 1 contains the values to count, column 2 is what the sequential field result should look like.

A             1

A             2

B             1

B             2

B             3

B             4

C             1

D             1

D             2

D             3

My plan was to approach this with an update cursor and the classic autoIncrement code, but it calculated all rows in the table with a sequential ID. I know the if/else is the wrong approach, but I was hoping the autoIncrement would start over within each if statement. I'm guessing I need to get a count of the unique values in Field1, then somehow calculate the sequential ID on the row as it counts? I saw a collections.Counter example, but wasn't sure if that could be applied here.

table = table path
fields = ('FIELD1,FIELD2')

rec=0
def autoIncrement():
    global rec
    pStart = 1
    pInterval = 1
    if (rec == 0):
        rec = pStart
    else:
        rec = rec + pInterval
        return rec

with arcpy.da.UpdateCursor(table, fields) as uc:
    for row in uc:
        if row[0] == 'A':
            row[1] = str(autoIncrement())
        etc.


texts: Why is there a difference in length for two lines created from the same points, one with ArcGIS Pro and the other with Geopandas?

texts:
Update
Question off-topic because I launched the ArcGIS treatment again and for an unknown reason this time I obtained the 21.56km length I was expecting.
It means the error comes from either my ArcGIS model or whatever the software did while executing.
In short :
What's the difference beween ArcGIS Pro's Points To Line and its Geopandas equivalent (using groupby() apply LineString) ?
On some data I get a different value of the line length, while on some other data the difference is insignificant.
Details on my case :
I have a pandas DataFrame containing about a thousand points with longitude and latitude, corresponding to a car drive recording.
I separated them in five zones (actually seven : from 0 to 6, but zones 0 and 3 are empty).
I want to get the length in kilometers for each zone.
Screenshot of the head of the table corresponding to zone 5 (I added the geometry column with the code, see next) :

What I used to do was using ArcGIS Pro :

Points to Line
then Add Geometry Attributes on the created lines
I would select "Geodesic Length" and the unit (km).

I'm trying to do the same on Geopandas.
Using the instructions from Turning GeoDataFrame of x,y coordinates into Linestrings using GROUPBY?, I got the part "transforming points into lines" covered.
Here is my code (still for zone 5) :
from shapely.geometry import Point, LineString, shape
import geopandas as gpd

geometry = [Point(xy) for xy in zip(df5.X, df5.Y)]
geo_df = gpd.GeoDataFrame(df5, geometry=geometry)

geo_df2 = geo_df.groupby(['Immatriculation'])['geometry'].apply(lambda x:LineString(x.tolist()))
geo_df2 = gpd.GeoDataFrame(geo_df2, geometry='geometry')

Then if I want the length in kilometers, I need to set a coordinate reference system and then change it to a projected coordinate system :
geo_df2.crs={'init' :'epsg:4326'}
geo_df2=geo_df2.to_crs({'init': 'epsg:3947'})

distance = [geo_df2.length[0]/1000]

I'm starting with the WGS 84 (EPSG=4326) because that's what ArcGIS uses for my data (both the points and line).
I chose to project to EPSG=3947 because it's what gives me the closest results to ArcGIS's results so far.
I realized that I got about 2 km of difference for zone 5 between ArcGIS's geodesic length and the one I get with Geopandas - given the precision I need, this is significative and I cannot ignore this error.

And we can see that it actually starts before the CRS projection :

Where does this significant difference for ONE zone and not the rest come from and can it be overcome?
It's the very same data for zone 5 in both cases so it doesn't come from the points.


texts: Gdal translate HDF5 to GeoTIFF library funtion

texts:
Im am trying to convert a radar file provided in HDF5 to a GeoTIFF using gdal.Translate() in my Python script.
An example of the HDF5 file can be downloaded here
I have already managed to get it working using the commandline and gdal_translate which gives me the expected result:
import os
os.system('gdal_translate -of GTiff \
-a_srs "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 y_0=0" \
-a_ullr 0.0 -3649999.1798970395 700000.9100000001 -4415003.157897039 \
"HDF5:\"RAD_NL25_PCP_NA_201910281110.h5\"://image1/image_data" \
"test-translate.tif"')

However, when I try to do the same using library function gdal.Translate() inside my script I get an error.
My code:
from osgeo import gdal

kwargs = {
    'format': 'GTiff',
    'projWinSRS' : '+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 y_0=0',
    'projWin' : [0.0, -3649999.1798970395, 700000.9100000001, -4415003.157897039]
}

fn = 'HDF5:\"RAD_NL25_PCP_NA_202101080825.h5\"://image1/image_data'
dst_fn = 'output.tif'

ds = gdal.Translate(dst_fn, fn, **kwargs)
ds = None

Gives me the following errors:
ERROR 0: -projwin_srs ignored since the dataset has no projection.
ERROR 1: Error: Computed -srcwin 0 -3.65e+06 700001 -765004 has negative width and/or height.

Am I doing something wrong?


texts: What is the formula for converting from WGS84 to EPSG:3071 (Wisconsin Transverse Mercator)

texts:
I have GIS software that can do this conversion but I need to know the actual formula / math involved. I'm writing code and need to perform a single conversion (WGS to EPSG:3071) so using a library, open source or not, involves too much overhead (20+ MB for one library). I don't need libraries of coordinate systems that will never be used. It would be best if I just write a function to perform the calculation.

How do I convert from WGS84 to EPSG:3071 (Wisconsin Transverse Mercator)?


texts: ArcGIS Desktop 10.5 - Create Line of Sight not showing intermediate points

texts:
I'm new to ArcGIS and in a college course.  I am trying to perform a Line of Sight analysis on a DEM raster, but it's only showing the endpoints, not the intermediate line with the visibility information in red and green.

I'm running it from the 3D analyst toolbar.

What it looks like

What the textbook says it should look like

I have double checked everything that the textbook, Discovering GIS and ArcGIS, has specified, and my prof has no clue what the error is.

What do you think is the problem here?

I know it's not much to go on, but I can't find a similar issue anywhere.


texts: How to use bbox in long/lat with google/bing/etc... using mapserver

texts:
The bbox of OpenLayers.Control.GetFeature works in utm when the baselayer is 900913 ie. google/bing/etc) and in long/lat when the baselayer is 4326 (ie. my postgis layers served by mapserver), no matter what overlays them. This is what openlayers says.

So, how can I force bbox to be in long/lat when I use one of my layers over google for example? I do want the bbox being in long/lat rather than 900913 (ie. using baselayer projection). I've read many webs about this but with geoserver, some posts about mapserver say that both EPSG:900913 and EPSG:4326 must be set in METADATA, I changed that and also assigned "srs" and "projection" (with the appropriate expression) in my layers (and in protocols) but with no effect in bbox at all. What else could I try?

Any idea will be appreciated, thanks in advance,

Gery


texts:
I am new in PostGIS and I am trying to convert a geometry field to GeoJSON Lat/Lng for visualization in a map as a layer.

I had a .HSP file that I exported to PostGIS with OpenJump program.

The .prj of .shp is:

PROJCS["ETRS_1989_UTM_Zone_30N",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-3.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

For get the SRID I searched in:
https://developers.arcgis.com/javascript/3/jshelp/gcs.html and I got the 4258 code.

But when I execute:

ST_AsGeoJSON(ST_Transform(ST_SetSRID(geo,4258),4326))

I get:

"{"type":"Polygon","coordinates":[[[526196.6063,4744503.4952],[526155.287,4744194.9938],[526067.3464,4744213.9994],[526088.0997,4744354.9984],[526109.542,4744513.4959],[526196.6063,4744503.4952]]]}"...

Which is not a Lat/Lng correct value.

What could be wrong?


texts:
Assume we have a road layer and we have a specific point and want to get all roads that surround this point or more precisely form the closest polygon around that point.
I miss a good approach - is there one?
(Also returning junctions closer to the point would be nice but optional)


texts: How to output a TIF with only 3 samples/pixel

texts:
I've been able to create a Topo from Canada's CanVec data. But after I 'Export as Image' and save my map to a TIF, I try importing it into Memory Map 6 and it says that it can't import the file because it has more than 3 samples/pixel.

After reading this page I wonder if this could be a color type issue (RGB or PCT)? Is there to change how many samples/pixel there is?


texts:
Very basic question, just started to use QGIS. I've been trying to import my own csv file with lat/long into QGIS. No luck, thought it was my excel file that was not with the proper format.

I tried with the file 'us_states_locations-1.csv' provided with QGIS when you download it. No luck either, same result. Here is the print screen. When I click ok, the system crash. Any idea of what I do wrong?


texts:
I am trying to convert a latlong point to UTM.

To define UTM projection, I need to calculate zone for the point.

I am trying to figure out the best way to do this.

One way to do this would be use longitude values to find the proper zone.

This would require lot of coding.

I am curious if there is a better way of doing this?


texts: Converting MapInfo TAB and/or MIF to SHP using free tools or ArcGIS for Desktop core product?

texts:
How do I convert MapInfo TAB and/or MIF files to SHP using only free tools or ArcGIS itself?

I've not had much luck trying to do this in the past and I don't want to have to shell out for the ArcGIS Data Interoperability extension which seems to be the 'official' way to go about this.


texts: Permanent Histogram Stretch with QGIS or FWTools

texts:
Is there a way to permanently stretch a histogram of a raster within QGIS or gdal?  I would like to stretch a satellite image with a Standard Deviation 2.5 stretch.  I would prefer to use QGIS or FWTools, but I am open to other open source tools.

Thanks...


texts:
I'm working in France.
I thought the TanDEM-X was a free and open DEM but that does not seem to be the case.

According to a few websites, the most famous free DEMs seem to be ASTER GDEM, SRTM, X-SAR and EU-DEM.

Their spatial resolutions appear to be equivalent to between 25 m and 30 m but what other parameters such as spatial accuracy, quality of modeling, etc define their usefulness?


texts: Displaying geopandas points with cartopy

texts:
I am trying to display a GeoDataFrame containing points, using cartopy.
I first tried to follow the indications given in the documentation of geopandas.
However, the points are not displayed:

Interestingly, if I use the Pandas internal plotting method, the points are indeed displayed, but the map extent is modified:

Here is the code:

import geopandas as gpd
import cartopy as ctpy
import matplotlib.pyplot as plt

map_data = '../data/stations.gpkg'

def plot_gpd(stations):
    # Reproj
    crs = ctpy.crs.Mollweide()
    crs_proj4 = crs.proj4_init
    geom_reproj = stations.geometry.to_crs(crs_proj4)
    # Plot
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection=crs)
    ax.add_feature(ctpy.feature.COASTLINE)
    geom_reproj.plot(ax=ax, markersize=3)
    plt.savefig('plot_gpd.png')

def plot_ctpy(stations):
    # Reproj
    crs = ctpy.crs.Mollweide()
    crs_proj4 = crs.proj4_init
    geom_reproj = stations.geometry.to_crs(crs_proj4)
    # Plot
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection=crs)
    ax.add_feature(ctpy.feature.COASTLINE)
    ax.add_geometries(geom_reproj, crs=crs)
    plt.savefig('plot_ctpy.png')

stations = gpd.read_file(map_data)
plot_gpd(stations)
plot_ctpy(stations)

Is there something I am doing wrong? Why is ax.add_geometries not working?

This related question is unfortunately not answering my problem: Add geopandas to cartopy


texts:
I've got a grid/raster in R that is represented as centroids of the grid instead of polygons. I am trying to merge that with points.

library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
grid_10000 <- st_make_grid(nc, cellsize = c(10000, 10000)) %>%
  st_sf(grid_id = 1:length(.),grid_prop = sample(c("one","two"),length(.),replace=T))
grid_centroids <- st_centroid(grid_50) %>% cbind(st_coordinates(.))
pts <- st_sample(nc, size = 5) %>% st_sf
# I can join the points easily like this
joined <- st_join(pts,grid_10000)

Now without grid_10000, how could I see which point gets which parameter simply based on the centroids. A method that could make a grid from the centroids would be best. The distances are unknown too. And 'sf' solution would be preferred.


texts: OpenLayers: How to move smaller vectors to front?

texts:
I've got a problem that's been bugging me for quite some time now. As you can see in the image below I've got som vectors that ends up beneath bigger ones, and despite all my efforts of trying to figure out how to use zIndex to solve this, I just can't seem to get it right.

What I want is, due to limitation in getting the data, to detect features within the area of the bigger ones, and then just change the zIndex or something in order to get the smaller ones "above" the bigger polygons, thus making them selectable.

Right now I'm pretty much at a dead end, and desperatly searching for any clue or solution to this.
Any ideas?


texts:
I tried to use the Terrain Profile tool, but it wont display any graph, no plot, just grey. It calculates the values in the table, but it seems some plugin is missing or not working but I have no idea what it could be.

Can anyone help me finding out the problem?

I am using Max OSX Lion

I installed matplotlib 1.3.0-3. several times


texts: Updating numerical column with autonumbers in sorted query table in ArcGIS for Desktop?

texts:
I am trying to find a suitable code that lets me to update a numerical column (namely, "RevOdrNo") of a sorted query table (not a separate table) with autonumbers (e.g. 1,2,3,4,5,6) in ArcGIS 10.3's Field Calculator without using the physical IDs (=not sorted).  I sorted this query table by right clicking the heading of the other numerical column (namely, "AX_Co") and clicking "Sort Ascending (or Descending)".

I have had a look at a number of the websites including the following link, but none of them seems to work for the sorted query table because they rely on the 'physical' order of the table rows.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00170000004s000000.htm (See the section of "Accumulative and sequential calculations")

I have had a thought about using the Sort tool in ArcGIS, but it asks for the output table name, which means that the query table will become a separate table, which I don't want to see.

Does anyone have any ideas of how to do it?


texts: ogr2ogr issues in clipping tab to shp where tab has points and lines together

texts:
I am running a simple script to clip and convert all tab files in a directory via osgeo shell

for %f in (*.TAB) do ogr2ogr "%~dpnf_Clip.shp" -clipsrc "clipped\AOI.TAB" "%f" -f "ESRI Shapefile"

It mostly works but for some files I get

  ERROR 1: Attempt to write non-point (LINESTRING) geometry to point
  shapefile.
  ERROR 1: Unable to write feature 602 from layer
  Water_Notes_DSC.
  ERROR 1: Terminating translation prematurely after
  failed

How can I get it to split these features into _point, _line and _poly automatically?


texts: I am converting data from spatialite to fileGDB (and reverse) for a project using QGIS and ArcGIS.  When I convert using ogr2ogr (or QGIS) from spatialite to fileGDB, the points are translated slightly during the conversion.  If I convert from spatialite to shapefile, there is no translation.  This is verified by doing a select by location query or by identifying the coordinates.  I have verified this condition in windows, and built gdal/ogr in linux using ESRI's fileGDB library.  Both produce the same results.  Based on my testing, it appears to be the write portion of ESRI's fileGDB library.  I have modified the XYSCALE and XYTOLERANCE settings, but neither appear to change the output.  Can anyone verify, or is there a solution?  The points are only translated slightly, typically 0.0001 ft, but it is enough to make spatial queries fail and prevent point updates.  GDAL version is 1.11.1.


texts: What is the compression gained by caching mod_tile's metatile instead png tiles?

texts:
What is the level of compression gained caching mod_tile's larger metatile instead of 256 by 256 .png tiles?

I am trying to determine the extra cost of storage resources for pre-rendering some set of .png tiles, which will be served as static files.


texts: Importing data from Excel into ArcGIS/QGIS

texts:
I have downloaded the hourly solar radiation data from NREL which is in CSV format but the problem is that the file does not have an exclusive column called latitude and longitude. There are 100s of CSV files like this and each file represents a latitude/longitude (we can say, point location). I have attached a screenshot of the file along with this post. I am not able to figure out how to import it in ArcGIS/QGIS.

The screenshot is for Bangalore city, India. The highlighted portion in red shows the latitude and longitude for the city. Can somebody please help me  import this in ArcGIS/QGIS?

Screenshot of the sample file (click to magnify):

Please find the same file attached here: http://www.sendspace.com/file/m9t19n

This file represents hourly solar radiation data for a particular latitude longitude (in this case, Bangalore).

Clarification: Essentially I would like to perform a site suitability analysis for solar for the state of Karnataka, India. My understanding is that, I will have point locations of the station representing the excel file in question, as the attribute. Then I will extract those locations with a certain DNI value, then intersect it with the land use to determine suitability.

I want to have the location of the station as well as all the attributes mapped to the location in the excel. So one location will have the entire excel as the attribute table (.dbf for the shapefile) and so on.