texts: How to create sequential numbers in a field by defining automatically the start, end and interval numbers?

texts:
I would like to create sequential numbers in a field called "virtual _Li".

I have used the field calculator using def autoincrement() function and manually calculate the interval number.
https://support.esri.com/en/technical-article/000011137

But I would like to have it in python.

The start point should be always the number in “Z_2018” that has “0” in “First_Dist” (Start = 1.006).
End number is the number in “Z_2018” that has the Maximum number in “FIRST_DIST” (End = 1.04).

The interval is:( (End) – (Start) ) / ( max in “FIRST_DIST” field/a constant number )

Constant number is always = 0.15

Start and end number can be + or – number.

In this case interval is 0.002875793

( (1.040) –  (1.006) ) / (1.773424/0.15) = 0.002875793

First create a field:

arcpy.AddField_management(ProfileGraph1WR1, "Virtual_Li", "DOUBLE", 10, 3)

then I need to create a “codeblock” and “Expression” which I don’t know How.


texts:
I'm using QGIS version 1.8.0-Lisboa and I'm trying to convert a vector line file and it's not working. I right click on the layer>Save as>Autocad file>DXF. After that I try to open the output file in Microstation (Bentley Software) and the file does not open at all. Can you please help me?

EDIT:
I already succeeded in converting to DXF, now I have another problem, I need a 3D file but the one I converted is just 2D, no height info

Can you help me get 3d dxf output?


texts:
I want to search for different buildings by using the labels I've added to them.

Like for instance, I want to search for the building with the label "1630", as in the picture. Is there any specific function or a plug-in for searching?


texts: Error inside pyinstaller app with pyqgis & pyqt4 on Windows

texts:
I am trying to generate a standalone exe using the pyqgis and related libraries.

I have borrowed from the app in the book 'Building Map Applications with QGIS', by Erik Westra and published by PACKT.

The code has been adapted in relation to this answer: Working Windows Standalone .EXE From QGIS Python Script Using PyInstaller?

Also the from xxxxx import * lines have been replaced and code updated accordingly. I have read from multiple sources that this can bother pyinstaller.

I have put all the code and associated bat files into a repository if anyone wants to clone it to try it out:
https://github.com/benboughton1/offlinemap

run.bat which is just running the app using standard python environment works fine.

pyinstaller runs OK too, but when I try to run the app I get the error:

Qt: Untested Windows version 10.0 detected!
LOADER: Running externalApp.py
Traceback (most recent call last):
  File "externalApp.py", line 4, in <module>
  File "c:\osgeo4~1\apps\python27\lib\site-packages\PyInstaller-3.2.1-py2.7.egg\PyInstaller\loader\pyimod03_importers.py", line 389, in load_module
    exec(bytecode, module.__dict__)
  File "qgis\__init__.py", line 26, in <module>
  File "c:\osgeo4~1\apps\python27\lib\site-packages\PyInstaller-3.2.1-py2.7.egg\PyInstaller\loader\pyimod03_importers.py", line 389, in load_module
    exec(bytecode, module.__dict__)
  File "qgis\PyQt\QtCore.py", line 28, in <module>
ValueError: API 'QDate' has already been set to version 1

Also the warning txt file has several missing modules:

missing module named PySide - imported by C:\OSGEO4~1\apps\Python27\lib\site-packages\pyinstaller-3.2.1-py2.7.egg\PyInstaller\loader\rthooks\pyi_rth_qt4plugins.py
missing module named org - imported by copy
missing module named qgis.core.QgsFeature - imported by qgis.core, qgis
missing module named qgis.core.QgsGeometry - imported by qgis.core, qgis
missing module named fcntl - imported by tempfile, subprocess
missing module named readline - imported by cmd, pdb
missing module named pwd - imported by posixpath
missing module named 'org.python' - imported by pickle
missing module named riscosenviron - imported by os
missing module named riscospath - imported by os
missing module named riscos - imported by os
missing module named ce - imported by os
missing module named _emx_link - imported by os
missing module named os2 - imported by os
missing module named posix - imported by os
missing module named resource - imported by posix


texts: Filling a field in a layer using a spatial relation with another layer/shapefile automaticly?

texts:
My question was kind of answered by @mgri at Filling a field in a layer using a spatial relation with another layer/shapefile?

If you have read that part, then I no longer need to repeat myself on this one. If not I will keep it simple by illustrate:

Just to note: I know how to create the process by using the following steps:

Create a buffer layer of my "bomen" point shapefile with "Fixed Buffer distance" tool.
Join this buffer layer to the straten shapefile and saving the records that intersect. Creating a new layer/shapefile in the process
Join the new merged layer/shapefile to the "Bomen" bomenschapefile and saving that overwriting the old shapefile.

How do I automate this so that when a new object is created in the "bomen" shapefile, this process gets executed?


texts: Adding a Layer from URL to ArcGIS Online with an API key in header

texts:
I am working with a vendor to get some realtime locations that they have posted as an API through their product added to an ArcGIS Online layer.

The API requires you to pass a key through after the call is made in a header. I have tried to pass in additional parameters

When I click add file from URL there is an option to add in parameters and I have attempted to pass in an "Authorization" with the API key, but it always shoots back an error.

Is it possible to add in an API key this way through ArcGIS online without using the ArcGIS API for JavaScript?


texts:
I created in Geoserver a store by using a PostGIS database (PostgreSQL 9.6) and I have some tables with geospatial data. When I add a new layer I cannot see the link 'Create new feature type' anymore. What am I missing?

I use Windows 10 x64, Geoserver 2.11.0, PostgreSQL 9.6.2 by BigSQL, PostGIS 2.3.


texts: PostGIS doesn't use spatial index with ST_Intersects

texts:
I imported the data from openstreetmap using osm2pgsql and copied it another table structure. If I made a query for containing rows using st_intersects, st_contain or someone else, postgis don't use the gist-INDEX. Always a seq-scan will be made. If I filter the rows first using the &&-operator (bbox of geometries), postgis uses the gist-index for the postgis functions. I don't understand why, where is the problem and how can I change this.

Here I've created a simplified table of the original (you can also use a planet_osm_polygon for queries):

CREATE TABLE osm_addr2 (osm_id bigint, geometry geometry);
CREATE INDEX osm_addr2_geom ON osm_addr2 USING gist (geometry);

Here I created a very simple example query:

EXPLAIN ANALYSE SELECT * FROM osm_addr2 AS addr
WHERE -- addr.geometry && (SELECT geometry FROM osm_addr2 WHERE osm_id=-332537) AND
  st_intersects(addr.geometry, (SELECT geometry FROM osm_addr2 WHERE osm_id=-332537));

Using st_intersects only, I have the following result (first line only):

Seq Scan on osm_addr2 addr  (cost=24.50..336.10 rows=387 width=40) (actual time=0.001..0.001 rows=0 loops=1)

If I uncomment the &&-operator the index will be used  (first line only):

Index Scan using osm_addr2_geom on osm_addr2 addr  (cost=49.15..57.41 rows=1 width=40) (actual time=0.031..0.031 rows=0 loops=1)

I'm using postgresql 9.3 with postgis 2.1, see postgis_full_version():

POSTGIS="2.1.4 r12966" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009" GDAL="GDAL 1.9.0, released 2011/12/29" LIBXML="2.8.0" LIBJSON="UNKNOWN" RASTER

I can reproduce this on Ubuntu 14.04 and Debian 7. I hope anyone can help because using the intersection operator will inflate my queries.


texts: Resurrecting contents of ArcGIS 9.x "samples" toolbox on 10.x?

texts:
Any ideas how to resurrect or uncover the contents of the ArcGIS Samples toolbox?

In version 9.x there were occasionally useful developer sample tools like Write Features To Text File and Batch Define Coordinate System among others. In version 10.0 the samples.tbx was deprecated, but still installed so that older scripts and models using the samples remain functional.

The .tbx is still there in 10.3, at ...\Desktop10.3\ArcToolbox\Toolboxes\Samples.tbx but marked as a Hidden File, and it doesn't show up in the toolbox list. If the hidden attribute is removed, it shows up. However looking at it in ArcCatalog/Map shows an empty container, yet the file on disk is 927kb (most of the tbx files are less than 100k).

I've verified on a 10.3 machine that the samples do indeed still work, so they're there. How can I see and inspect them directly?

I need to write arcinfo ungenerate files, and there's a bug in either the code or the docs for WriteFeaturesToTextFile, I want to fix it. While I'm there I might as well let some other possibly useful code see the light of day again.


texts:
I am new to programming in ArcMap, and am using ArcMap 10.5.
When I use the tool Select (Extract-Analysis Tool), the result is what I expect.
However, When I run the Phyton Snipped code ArcMap displays:
“Parsing error SyntaxError: invalid syntax (line 1)”
I have searching the web with no success. Can you tell me how to fix this error?
On the other hand, I have to apply the tool 130 times to export the polygons shapes one by one. How can I include the code in a for cycle?

arcpy.Select_analysis(in_features="D:/AKKGIS/cogo_plgs/001_130_cg_plgs.shp", out_feature_class="D:/ AKKGIS/img_clips/001_clp.shp", where_clause='"ppg_id" = '001'')


texts: Creating NDVI time series from single polygon of shapefile using Google Earth Engine

texts:
I am new in Google Earth Engine, I am trying to plot ndvi timeseries from a single polygon of shapefile, but when I run the code I get "Error generating chart: Collection query aborted after accumulating over 5000 elements".

var broad_leaved= ee.FeatureCollection('users/spatola/MCD45A1_broad_leaved_311/MCD45A1_2004_clc311');
var modis = ee.ImageCollection('MODIS/006/MOD13Q1');
var modisNDVI = modis.filterBounds(broad_leaved)
                    .filterDate('2002-01-01', '2009-12-31')
                    .select('NDVI');
                    print(modisNDVI);
var timeseries = ui.Chart.image.seriesByRegion(
  modisNDVI, broad_leaved , ee.Reducer.mean(), 'NDVI', 250)
  .setOptions({
    title : 'vegetation series',
    xAxis : {title : 'YEAR'},
    yAxis : {title : 'NDVI'}
    });
  print(timeseries);


texts: Load 3D shapefiles in QGIS 3.0 where Ground Wall Roof Solid are separated

texts:
The brussels region distribute 3D shapefiles in opendata.

http://urbisdownload.gis.irisnet.be/fr/download/

The spatial objects are separated in many files :

Ground
Wall
Roof
Solid

There is no 3D informations about 3D in attributes.

AND the geometries are flat:

How to render 3d?


texts:
I'm using QGIS 3.10.0 A Coruna.
As you see in the picture, I have an unregular grid with some buildings on them. I would like to fill out all the other rectangles between the points with a polygon automatically. That means every 4 points should create a seperate polygon. Is there a possibility to do so?


texts: ArcGIS Server 10.1 and Postgres Geoprocessing Service issue

texts:
I'm hitting an error when publishing a model from ArcGIS Desktop 10.1 to ArcGIS Server 10.1.  The model I want to publish uses a layer from a Postgres DB defined via a Query Layer in my MXD.  The model is very simple for now, run the Select function on the Query Layer and return the result as a shapefile.

Publishing fails with a "Failed to create service definition" message and a "ERROR 001270: Consolidating the data failed" message in the results window.

This error links to http://resources.arcgis.com/en/help/main/10.1/index.html#//00vp0000002r001270

I have taken all of the steps mentioned therein - shortening service names, paths etc, and ensuring that the DB connection is registered in my ArcGIS Server instance.  I can add the PG layer as a standalone map service which works fine.

Thanks for any help or suggestions!

EDIT:
I should add that my Postgres instance (9.1.3) is not configured as an ESRI Enterprise Geodatabase, so that could be an issue?  PostGIS 1.5.3 is installed on the DB, which is currently used for web mapping (MapServer) and hopefully data connections into ESRI GP services.


texts: While logging into GeoServer the username and password are not taken the ones which I had given during installation. Any suggestions?


texts: Use python standalone script as QGIS tool

texts:
I have python scripts that I use to do some geo-processing. I need to share it with QGIS users that are not familiar with python. I've written QGIS tools in the past using the python interface from QGIS, but the scripts I want to use now contain some modules (e.g., pandas) that I was unable to install on QGIS.

Is there a way to convert my python scripts "as is" to run in QGIS as a tool for non-python users? I can also convert my script to .exe if it makes it easier.

Edited: My inputs are vector file as geojson and my rasters are geo-tifs. I was able to install some modules in QGIS in the past but could not install others. I'm looking for a straight-forward solution (if it exists). Meaning, I have a python (standalone) script that works great on IDE and also works great as exe file (from the cmd), I want to take this script "as is" and convert it so that it would work on QGIS. Is there an easy and simple way to do that?


texts:
I just looked for MapTips in ArcGIS Pro 1.2 and have not been able to find them in either the application or its documentation.

Has anyone found them in ArcGIS Pro 1.2, or know when/whether they are expected to be added?

For those unfamiliar with map tips, the ArcMap documentation describes them:

  MapTips pop up as you hover the pointer on a feature, providing a
  quick way to see the name of that feature or a particular piece of
  information without having to use the Identify tool or HTML pop-up
  windows. ArcMap lets you choose which attribute field you want to
  display as your MapTip.

Personally, I do not think that they are overly important to have, and am happy to just label or identify/popup to get the same/more info about visible features.


texts: Is it possible to access ArcObjects with IronPython?


texts:
I'm using ArcGIS 9.3.1 and trying to get the code to create a new feature class in an existing personal geodabase to work but I am having no luck. When running in debug mode from within vb.net it will fail with the standard "LoaderLock" error where it tells you to not run managed code inside a DllMain. I'm sure this is a simple problem that I'm just not seeing, but I've been searching the forums for a couple of days not with no luck. If anyone has insight after looking at the code below I would be very thankful.

Dan

Dim pGeomDef As IGeometryDef = New GeometryDef

Dim pGeomDefEdit As IGeometryDefEdit = CType(pGeomDef, IGeometryDefEdit)

With pGeomDefEdit

.GeometryType_2 = esriGeometryType.esriGeometryPolygon

Dim pSpatRefFact As ISpatialReferenceFactory2 = New ESRI.ArcGIS.Geometry.SpatialReferenceEnvironment

Dim pGeoCoordSys As IGeographicCoordinateSystem = pSpatRefFact.CreateGeographicCoordinateSystem(ESRI .ArcGIS.Geometry.esriSRGeoCSType.esriSRGeoCS_WGS19 84)

.SpatialReference_2 = pGeoCoordSys

End With

Dim pFields As IFields = New Fields

Dim pFieldsEdit As IFieldsEdit

Dim pField As IField

Dim pfieldEdit As IFieldEdit

pFieldsEdit = pFields

pField = New Field

pfieldEdit = pField

With pfieldEdit

.Name_2 = "ObjectID"

.Type_2 = esriFieldType.esriFieldTypeOID

End With

pFieldsEdit.AddField(pField)

pFieldsEdit = pFields

pField = New Field

pfieldEdit = pField

With pfieldEdit

.Name_2 = "SHAPE"

.Type_2 = esriFieldType.esriFieldTypeGeometry

.GeometryDef_2 = pGeomDef

End With

pFieldsEdit.AddField(pField)

Dim pWorkSpaceFactory As IWorkspaceFactory = New ESRI.ArcGIS.DataSourcesGDB.AccessWorkspaceFactory

Dim pFWS As IFeatureWorkspace = pWorkSpaceFactory.OpenFromFile("c:\Dan\Data\DataFi tness-1AOI.mdb", 0)

Dim pFeatClass As IFeatureClass

Dim pUID As New ESRI.ArcGIS.esriSystem.UID

pUID.Value = "esriGeoDatabase.Feature"

'FAILING HERE

pFeatClass = pFWS.CreateFeatureClass("TestFC1", pFieldsEdit, pUID, Nothing, esriFeatureType.esriFTSimple, "SHAPE", "")


texts:
The loop runs once, then i get this error message.

Runtime error  Traceback (most recent call last):   File "<string>", line 17, in <module>   File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\management.py", line 6997, in MakeFeatureLayer     raise e ExecuteError: ERROR 000725: Output Layer: Dataset out_VHQ_7 already exists.

To make this code work i think i need to add something to the outputname "out_VHQ_7" to avoid that when the loop runs the second time it stops because it cannot overwrite the output that has been createt in the first run of the loop

>>> FieldToUse = "Symbol" # change this to your field name with the unique values
... for fc in fcList:
...     uVals = [] # new empty list
...     # loop through each feature with a searchcursor
...     with arcpy.da.SearchCursor(fc,FieldToUse) as SCur:
...         for row in SCur:
...             # check if the value is already in the list
...             if row[0] not in uVals:
...                 uVals.append(row[0]) # add new value to list
...
...     #loop through each unique field value
...     for FieldValue in uVals:
...         # create a layer with a where clause, this is for numeric fields
...         arcpy.MakeFeatureLayer_management(fc,"out_VHQ_7","{} = \'{}\'".format(FieldToUse,FieldValue))
...         # use "{} = \'{}\'".format(FieldToUse,FieldValue) to quote the value
...         # if the input field is a string field
...         OutTIFFname = "{}_{}_kd.tif".format(os.path.splitext(fc)[0],FieldValue)
...         outKDens = KernelDensity("out_VHQ_7","NONE",4670)
...         outKDens.save(os.path.join(outworkspace,OutTIFFname))


texts:
I am working with INSAT3D Images from mosdac for air quality monitoring as my final thesis. I had downloaded the data from mosdac which are in HDF format and i tried to convert it to GeoTiff using HEG tool, but the operation failed since the pixel size which have to give in  the geotiff for the conversion was incorrect.

Is there any software for the conversion of HDF image for INSAT3d to a Geotiff file other than HEGtool?

If not can anyone let me know the pixel size which have to give in  the geotiff for the conversion?


texts: Using ArcGIS, can data relationships be created between separate geodatabase containers?

texts: I have several ESRI geodatabases that contain specific detailed information around a central theme.  Without revealing the specifics about the theme, each geodatabase has its own domain-driven subtypes and value tables, but there appears to be no attribute connection between geodatabases.  Presuming a relationship - other than geospatial - can be surmised between databases, how would this relationship be created?  I am using ArcGIS 10.2.


texts:
I'd like to reconfigure an ArcGIS Server 10 Java implementation (single machine, Windows Server 2008 R2) and I plan to re-run the post-install.

Complete: I do not need to save anything on the current machine. To get the cleanest reinstall, I intend to delete the SOM, SOC & "Webservices" accounts; the agsusers & admin groups, and the five directories created (arcgiscache, arcgisindex, arcgisinput, arcgisjobs, & arcgisoutput).

Any other suggestions for redoing an ArcGIS Server 10 post-install?

Update: I reran the post-install and have quite a few notes. I posted them below and made this a community wiki.


texts: Failure of non existant SQL statement after executing statements from a file with ArcSDESQLExecute object

texts:
As you see in  my code I iterate and execute all the sql statements in an sql file using ArcSDESQLExecute object.

It does work but when all statements have been executed the object returns 'None' (and thus the elif sde_return == False: part) indicating a failure and I don't know why.

What am I doing wrong?

I've used the first example in Executing SQL using an ArcSDE connection to build my code.

My code:

import arcpy

try:

  sqlFileSt = arcpy.GetParameterAsText(0)
  destSDE = arcpy.GetParameterAsText(1)

  arcpy.env.workspace = destSDE
  sde_conn = arcpy.ArcSDESQLExecute(destSDE)

  fd = open(sqlFileSt, 'r')
  sqlFile = fd.read()
  fd.close()
  sqlFile = "".join([s for s in sqlFile.strip().splitlines(True) if s.strip("\r\n").strip()])
  sql_llist = sqlFile.split(";")
  for sql in sql_llist:
    try:
      sde_return = sde_conn.execute(sql)
    except Exception as err:
      arcpy.AddMessage(err)
      sde_return = False

    if sde_return == True:
      arcpy.AddMessage("SQL: {0} executed.".format(sql))
    elif sde_return == False:
      arcpy.AddMessage(arcpy.GetMessages(2))
      arcpy.AddMessage("SQL: {0} failed.".format(sql))

except Exception as err:
    print(err)
    arcpy.AddMessage(arcpy.GetMessages(2))

And the contents of in my sql file:

CREATE TABLE "APP_UPDATES"
  (
    "OBJECTID" NUMBER(38,0) NOT NULL ENABLE,
    "DATE_ACT" DATE,
    "NAME_TABLE" NVARCHAR2(50),
    "ELEMENTS_DEL" NUMBER(*,0),
    "ELEMENTS_ADD"   NUMBER(*,0),
    "RESULT" NVARCHAR2(15),
    "ERRORS" NCLOB
    );

CREATE TABLE "APP_CONFIG"
   (
    "C_KEY" VARCHAR2(50 BYTE) NOT NULL ENABLE,
    "C_VALUE" VARCHAR2(2048 BYTE) NOT NULL ENABLE,
    "COMMENTS" VARCHAR2(255 BYTE)
    );

CREATE TABLE "APP_QUERIES"
   (
    "NAME" VARCHAR2(255 BYTE) NOT NULL ENABLE,
    "QUERY" VARCHAR2(2048 BYTE) NOT NULL ENABLE,
    "NOTES" VARCHAR2(255 BYTE)
    );


texts: Using peewee ORM with spatialite databases?

texts:
I would like to use the peewee orm in front of some Spatialite databases.  Getting peewee to connect to a Spatialite database was relatively simple, and basically boiled down to:

class SpatialiteDatabase(SqliteExtDatabase):
    def __init__(self, *args, **kwargs):
        super(SpatialiteDatabase, self).__init__(*args, **kwargs)
        self.load_extension('mod_spatialite')

    def initialize_connection(self, conn):
        LOG.warn('initializing spatialite metadata')
        self.execute_sql('select InitSpatialMetadata(1);')

playhouse.db_url.schemes['spatialite'] = SpatialiteDatabase

So, great, now we can get a connection to a database with the Spatialite
extension enabled.  But I'm not clear how to create Peewee models with
geometry columns, since unlike in, say, PostGIS, Spatialite requires
the use of particular SQL functions (like AddGeometryColumn) to
manipulate tables.

Ultimately, I would like to be able to write a model something like this:

class Placemark(Model):
    id = PrimaryKeyField()
    name = TextField()
    alias = CharField(20)
    geom = GeometryField(srid=4326, geom_type='POINT')

And get the equivalent of:

CREATE TABLE "placemark" (
  "id" INTEGER NOT NULL PRIMARY KEY,
  "name" TEXT NOT NULL,
  "alias" VARCHAR(20) NOT NULL,
);
SELECT AddGeometryColumn("placemark", "geom",
    4326, "POINT", "XY");

My best option so far has been to populate the geometry fields manually after the fact.


texts: QGIS 3.18 and 3.16 LTR install. GRASS and SAGA not working. MacOS OSX

texts:
I have tried to update QGIS on Mac to either version 3.18 or 3.16 using the installers on the QGIS page (https://qgis.org/en/site/). I have been having issues getting the GRASS and SAGA providers running.
Error message below.
GRASS init error: Problem in GRASS initialization, GRASS provider and plugin will not work : Module built against version 2021-02-10T17:13:06+00:00 but trying to use version 2021-02-10T17:08:53+00:00. You need to rebuild GRASS GIS or untangle multiple installations.
No idea how to rebuild GRASS or untangle installations
Big sur 11.2


texts:
Can I add a photogallery in the Cartodb Infowindow?

I know how to add an image (single) in the infowindow, but can a photogallery be added as well? Small thumbnails that become larger when clicked?

Can this be done with an iframe or something like that?


texts: Connect a point geometry with a graphics label on a map

texts:
I am creating a map with traffic signals in a city.

Creating a point layer and adding my symbols as SVG mark is appropriate for 90%. But some signals explain alternativ routes used during events. These are not legible anymore when size is small enough to fit the map.

That's why I'd like to create a line between the point and the picture, similar as it is possible with labels. This would allow me to place these signals more outside and in a reasonable size.

Any ideas?


texts:
I am a beginner with PostGIS.

I would like to know if it’s possible to use Count or “>” with ST_Intersect.

Which districts sit in 3 neighbourhoods ?

SELECT district.name
FROM district
JOIN quartier ON district.quartier= quartier.district
WHERE count (ST_Intersects (district.geom, quartier.geom))


texts: QGIS 3x GDAL raster to vector (polygonize) option doesn't work

texts:
I want to convert my Geotiff map into the shapefile.
For this purpose, I used the GDAL->Raster Conversion -> Polygonize option.

I have the status "Processing algorithm" for nearly 2 hours!
I found the similar problem here:
Long Time for conversion raster file to vector in QGIS
but the guy had the file over 1GB. I already clipped my Geotiff map trimming down the size to 30MB.
Why does it take so long then?
Have I missed something?
Or shall I use other options...
Converting Raster to vector and creating polygons based on each pixel using QGIS


texts: PostGIS resampling two rasters to common grid system

texts:
I am testing the PostGIS map algebra (ST_MapAlgebra) on two different raster tables.

Each table contains a single image/row. They have the same SRID but different, overlapping extents and different ground resolution. What I did was:

SELECT ST_MapAlgebra(t1.rast,t2.rast
        ,'([rast1]+0.001)/([rast2]+0.001)') AS rast
FROM t1, t2

And I got an error saying that

  rt_raster_iterator: The set of rasters provided (custom extent
  included, if appropriate) do not have the same alignment

Now, what I want to do is to resample t1.rast and t2.rast to their common area (i.e. intersection) at the higher resolution of the two rasters.

Is there a function for this in PostGIS (2.1+)? and, if not, how should I go about making such a query?

I looked at ST_REsample, and it seems that it can only sample one image to the reference grid of another. Also, the documentation says that it does not handle SRID since PostGIS 2.1.


texts: Sum Values for multiple unique field values

texts:
I am trying to modify an existing table that I can relate to a feature in Arcmap(v10.5).

Currently, the table fields and values look like this:

I would like to summarize the table.  For each unique "property number" value, I would like to have the table list unique "fagtype" values.  And for each "fagtype value", I would like the table to sum the "fagacres" for each value.  I would like the "facres" and "fagirrig" fields to just list unique values as well.

I found a starting place under the "search cursor" summary on the ESRI help page.  But it's a very limited start and I'm not sure where to go from here, or even if a searchcursor is the best route.

import arcpy

agtable = arcpy.GetParameterAsText(0)

fields = ['PID', 'fagtype', 'fagacres', 'facres', 'fagirrig']

with arcpy.da.SearchCursor(agtable, fields) as cursor:

    for row in cursor:


texts:
I have written a python script which has one output line-vector layer. The python script should calculate some value and then the script should set the attribute (with that value) of a newFeature (which will be added to the output layer).

As I have tried it with code from QGIS2.18, I get the following error:

  Traceback (most recent call last):
  File "", line 125, in processAlgorithm
  KeyError: '0'

My code looks like this:

# loop through edges of inLayerEdges and merge
    for feat in inEdges.getFeatures(QgsFeatureRequest()):

        # create new edge(line) for outputLayer
        newEdge= QgsFeature()
        currentPriority = feat.attributes()[idxPriority]
        currentEdgeId = feat.attributes()[idxEdgeId]

        if currentEdgeId.find('From') == -1:
            edgeIdTrimmed = currentEdgeId.replace('To', '')
        else:
            edgeIdTrimmed = currentEdgeId.replace('From', '')

        if edgeIdTrimmed in listOfEdgeIds:
            # if current edge id already processed -> skip it.
            #print "current edgeid already in list. Continue ..."
            continue

        else:
            # if not, find all entries in the graph which share the same current edge id
            # and merge them
            listOfEdgeIds.append(edgeIdTrimmed)

            newEdge.setAttribute(EDGEID, edgeIdTrimmed)
            newEdge.setAttribute(PRIORITY, currentPriority)

    (...)

    # add the feature to the layer
    (res, outFeats) = outdp.addFeatures([newEdge])


texts: Eventually, I was able to get the GDAL installed and compiled correctly using Visual C++ 2010. I even tested some of its methods on some geo-referenced raster files, and it looks like that all is working fine. However, I am interested in running GDAL with some sort of hardware acceleration? Does anyone know if there is any support for GDAL with some acceleration libraries (e.g. CUDA/OepnCL, or maybe even with using OpenGL for the sake of GPGPU, or SIMD)?


texts:
In ArcGIS 10, I have several maps describing both habitat characteristics and climate. I also have a raster map of presence/absence for a species that I wish to examine using my other maps. All my maps are in raster format. I wish to merge all of the information from all of my individual maps into one table, that I can use for regression analysis, testing what features of the habitat and climate my species depend on.

I have first ensured that all of my maps have the exact same extent and shape (using the Clip tool).

I have also created a fishnet (my polygon), and wish to calculate the mean of each of my single maps within each of the grids of the fishnet (except for my species, where I wish to find the majority value for each cell). I used the Zonal statistics as table to extract the majority value of my species dataset into the fishnet. In order to check that the data was properly transferred to the table, I used the Join Field Tool to join my zonal table to the fishnet polygon.
To check, I used Select by attribute and choose all instances where my species is present (value 1). The problem is that the chosen part of the dataset does not correspond with the actual species distribution (see attached picture). The picture shows the selected part of my fishnet (blue) compared to the actual distribution (presence = black, absence =gray). If done correctly, I assume all black and blue dots should overlap.

I have also attached a picture of my model and one of the specifics in the Zonal statistics as table Tool.

I hope someone can tell me what I have to do in order to join all of my dataset into one table in a correct way.


texts: Making QGIS GUI responsive while adding multiple layers in Plugin

texts:
I have a tool button that add many layers from python code (WFS, PostgreSQL, ...), about 50 layers (30s to load) in QGIS.
During the process the application seem to be crashed "not responding". I've tried to add layers with a QgsTask class but it crashes too because i have to be able to access to QgsProject instance. In the documentation there is an avertissement about this, so i can't use a QgsTask object.
Is there an other solution to keep the GUI alive when I make a heavy work?
EDIT :
I have isolated the part of my code for adding layers
def _loadLayers(self,layers, grp) :
    for layer in layers:
        currentLayer = QgsVectorLayer(uri.uri(False), layer['name'], layer['type'])
        ...
        if (currentLayer is not None) :
           currentLayer.triggerRepaint()
           qgis.utils.iface.layerTreeView().refreshLayerSymbology(currentLayer.id())
           QgsProject.instance().addMapLayer(currentLayer, False)
           node = grp.addLayer(currentLayer)

conf  = Configuration.getOptions() # Json file with layers definition
if (conf is not None) :
   if ("groups" in conf) :
      ltr = QgsProject.instance().layerTreeRoot()
      for child in conf['groups'] :
          if ("groups" in child) :
              self._loadGroup(child, ltr)
          elif ("layers" in child) :
              if (ltr.findGroup(child['name']) is None) :
                 grp = ltr.addGroup(child['name'])
                 self._loadLayers(child['layers'], grp)


texts: locator service to use with ArcGIS Javascript API v2.1

texts:
I am hoping to find a locator service that will locate by address or intersection, and that will work with version 2.1 of the ArcGIS Javascript API.

I tried this locator given in the examples, which works well with 3.1, but it doesn't seem to work with 2.1 :

http://help.arcgis.com/en/webapi/jav...r_address.html

Unfortunately I think I need to stick with v2.1.  SO far I haven't had any luck.

Thanks in advance, Jason


texts: Generating polygon representing rough 100km circle around latitude/longitude point using Python?

texts:
Given a latitude/longitude point anywhere on earth, I would like to assemble a polygon that roughly resembles a circle with a 100km radius around that point.

It doesn't have to be pin-point accurate, just good to within a few hundred meters. I'd like to control the number of points returned (20 is probably fine).

I'm already using the Shapely library for Python, but it looks like that can't quite get the job done on its own.

My ideal function would look something like this:

polygon = circle_around_point(
    latitude=29.9500,
    longitude=-90.0667,
    radius_in_meters=1000,
    num_points=20
)

polygon would then be a list of 20 lat/lon points around the edge of that circle.


texts: Rotating non-symmetrical symbols, or avoiding upside down symbols

texts:
I have a lovely symbol of a vehicle on a map (as an SVG path). I'd like to show current heading.

Trouble is, cars aren't symmetrical. I am ending up with a lot of vehicles that look like this. Not what my bosses want to see.

My question: is there a way to flip this symbol or load a different symbol if the heading would make it be upside down? Or, more succinctly, how to I keep the cars right side up?

My code:

    tempRenderer = new UniqueValueRenderer(basicSymbol, "age");
    tempRenderer.addValue(0, new SimpleMarkerSymbol(SimpleMarkerSymbol.STYLE_PATH, 20, '', new Color([0,0,0,1])).setPath(truckSVG));
    tempRenderer.setColorInfo({
        field: "temp",
        stops: [
            { value: -30, color: new Color([255, 255, 255]) },
            { value: 0, color: new Color([200, 128, 200]) },
            { value: 30, color: new Color([128, 128, 255]) },
            { value: 60, color: new Color([128, 255, 255]) }
        ]
    });
    tempRenderer.setSizeInfo({
        field: "age",
        stops: [
            { value: 0, size: 40 },
            { value: 1, size: 10 },
            { value: 30, size: 4 }
        ]
    });
    tempRenderer.setRotationInfo({
        field: "heading",
        type: "geographic"
    });
    truckLayer = new GraphicsLayer();
    truckLayer.setRenderer(tempRenderer);
    map.addLayer(truckLayer);

As you can see, I'm showing the last 30 minutes of data for where the vehicle has been. The newest one gets a cool truck icon, while older ones get a dot. The temperature at the vehicle's location determines the color, the age of the data determines the size of the dot.

The data is being pulled in as a plain old json file (not an esri service). I suppose it'd be possible to put the latest point in a separate layer.

I can create an SVG facing the opposite direction if necessary and manipulate the headings and set a flag if needed.


texts: QGIS - Export raster layer to SVG with maximum detail

texts:
I'm trying to export the "Stamen Terrain labels" raster shapefile to an SVG.
So I added it to a new QGIS project then I'm trying to use print composer to export the SVG.
I want the entire world with lots of data when I zoom in.  So I zoom in to the layer in QGIS before opening the print composer to get the level of detail I want, but then when I go to the print composer, I'm only getting the portion of the map that is visible on screen.
If I scroll out until the entire layer is visible on the screen, then when I go to Print composer I don't get the level of detail I want for the labels.  For example, I'll just have the continent names, not the country names / big cities, etc.
So how do I get the level of detail I want but for the entire world loaded into the print composer so I can then export it to a SVG?  It seems like it should just be a "canvas size" setting someplace, but I don't see it.


texts: Is there a better geoprocessor for viewshed analysis in ArcGIS JavaScript API?

texts:
I am working on a viewshed viewer application for webcams. Viewshed variables that are critical to my application include: elevation, azimuth, and direction.
ArcGIS' own Viewshed Analysis example uses this geoprocessor. But the variables it takes seem to only include a point and a max viewing distance...
var gpUrl = "https://sampleserver6.arcgisonline.com/arcgis/rest/services/Elevation/ESRI_Elevation_World/GPServer/Viewshed";
var gp = new Geoprocessor(gpUrl);
        gp.outSpatialReference = {
          // autocasts as new SpatialReference()
          wkid: 102100
        };
var params = {
            Input_Observation_Point: featureSet,
            Viewshed_Distance: vsDistance
          };

          gp.execute(params).then(drawResultData);

Is there a way to input my desired variables in this existing geoprocessor? Or is there something even better in the API?


texts: Pasting point data into PostGIS with M Values

texts:
I have point data exported from Trimble Insphere as .shp. I am trying to paste the data into our PostGIS database but keep getting:

PostGIS error while adding features: ERROR:  Column has M dimension but geometry does not

I have tried exporting to PointZM (edit: export maintains native point 25d geometry), set M, and a few other things with no luck. I don't have admin privileges for the database so I need to work it out on my end if possible. Any ideas?

Edit: Using QGIS 3.6


texts: How to extract certain values from a table using the SearchCursor Tool

texts:
I am trying to write a python script to extract certain values from a table: the table I am referring to is a big collection of nitrate values for different water depths, which are deposited in the columns of the table. As I only need the value of the surface and the deepest point, I want to search through the rows and extract the last value that is not 0. I have started writing a script using the SearchCursor Tool but get stuck at the point, where I want it to search for the first 0-value and then go back and print the value fro mthe column before... Does anyone have an idea how to solve that problem?

import arcpy

# Set the Workspace
arcpy.env.workspace = "D:\Teresa\Kerstin\SouthernOcean\03_workspace\Teresa"

# Make table
table = "C:/Users/theidema/Desktop/OxzUti_GridP_Annual.csv"

#Create the search cursor
cursor = arcpy.SearchCursor(Table)

#Iterate through the rows
row = cursor.next()
while row:
    print (row.getValue(field))
    row = cursor.next()

Here is a Screenshot of the table (depths go down until 5500M)


texts: Finding data driven pages for ArcGIS Pro?

texts:
I am working on importing .mxd files into ArcGIS Pro 1.2.  Everything appears to import smoothly except I lose my data driven pages. I read Pro 1.2 is supposed to support them.  Yet, when I do a geoprocessing search in Pro, nothing specific to data driven pages appears, and it is not located under Cartography Tools like it is for Desktop.

So, my questions is really broken down into three parts.

Does Pro 1.2 support data driven pages?
If so, where is it located?
Is there a way to import data driven pages into Pro from ArcMap?


texts: Snapping multiple polylines to points nearby

texts:
I have a network of polylines and nodes. Each polyline starts at a point and ends at other. To be more specific, the network consists of links and nodes.

The problem is that, the polylines are slightly (maybe feets) offset to the points. Is there any way that I can snap all the ends of the polylines to the nearest point??

Now the new network would have points exactly at ends of polylines.


texts: How can I determine whether a feature class is a 3D feature class (HasZ)?

texts:
I would like to be able to determine whether a feature class passed to a procedure is a 3D feature class (hasZ) or not.
How can I do this in VB .Net & ArcObjects?


texts: Syntax Error related to File Path / Name in ArcPy Script

texts:
I'm using ArcGIS Desktop 10.2.2.

I'm getting the following error when I run a script in Arc.

SyntaxError: EOL while scanning string literal

The error message points to the line where I set my workspace variable "wspace".

Here are the relevant portions of my code:

# Set the local variables
wspace = r"Q:\GW\EC1210WQAEH_QESEA\CSSP_ATL\Users\vanvulpenc\CSV_to_Shapefile\temp_workspace"
x_coords = "Lng"
y_coords = "Lat"
output_Layer = "output_Layer"
saved_Layer = "To_Shapefile.lyr"

Then later in my script the one line that references the "wspace" variable is:

arcpy.CopyFeatures_management(wspace+saved_Layer, Output_FC)

Is there a better way to format my path perhaps?


texts: Setting Geometry Fail Error on ESRI server

texts:
I'm trying to insert a point on a feature layer on a ESRI server.
Sometimes it works nicely, but sometimes I get an error when I send the save instruction:
"Setting of Geometry Fail"

This happens when I build the geometry from a mappoint. When I compare the geometry from a mappoint with a normal geometry, I don't find any difference.
The mappoint is made only by using X,Y and WKID.

Can anyone have a clue of what I'm doing wrong?

Thanks


texts:
I have a polyline represented by LineString class. I need to add additional points in between its points at a distance, say d. So if P1 and P2 are first 2 points in LineString, I would like to have points inserted between P1 and P2 such that distance between points is d.
Does locationtech.jts offer such an API?


texts: Subsetting a shapefile and saving it using PyQgis

texts:
I am very new to PyQgis.

I have a shapefile containing a lot of information for each place, it looks like this

id    stat_1   stat_2 ... stat_n
i1     0        0          2
i2     1        1          0
            ...
iT     0        0          2

I want to create n shapefiles, containing for each of them: id and sta_i and only the locations which statistic is not equal to 0. Something like:

id    stat_i
i3     1
ij     2
            ...
iT     4

Doing this manually would be quite tedious, which is why I decided to get started with PyQgis.
I figured out how to select my main layer and how to select the features I want ( I actually just figured out how to print them).
How can I actually create a new shapefile with the subsetted locations and just the two desired fields?
I will sort out how to loop over the fields later.

layer=None
fields=["id","stat_i"]

for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "Layer_Name"
        layer = lyr
        break

request=QgsFeatureRequest().setFilterExpression(u'"stat_i"!=0')
it=layer.getFeatures(request)
for feat in it :
        print feat.id()


texts:
Here is my code for getting a WMS layer in OpenLayers from Geoserver.

var wfs2 = new OpenLayers.Layer.Vector("WFS", {
                    strategies: [new OpenLayers.Strategy.Fixed()],
                    protocol: new OpenLayers.Protocol.WFS({
                    version: "1.1.0",
                    url: "http://localhost:8080/geoserver/wfs",
                    featurePrefix: "topp",
                    featureType: "states",
                    featureNS: "www.openplans.org/topp",
                    srsName: "EPSG:4326",
                    geometryName: "position"
                  })
         });

So i get this problem in developer tools for chrome :

XMLHttpRequest cannot load http://localhost:8080/geoserver/wfs. Origin null is not allowed by Access-Control-Allow-Origin.

in developer tools for chrome

I saw that people recommend to set up a proxy, but my .html is in the geoserver/data_dir/www folder, and don't call external urls.

I however tried it, so I created the folder cgi-bin in /usr/local/geoserver/cgi-bin/, where i put the proxy code I found here :http://trac.osgeo.org/openlayers/wiki/FrequentlyAskedQuestions#ProxyHost

And then I added in the cript of my .html this : OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

But I then get this error in developer tools :

Failed to load resource file://localhost/cgi-bin/proxy.cgi?url=http%3A%2F%2Flocalhost%3A8080%2Fgeoserver%2Fwfs

Probably the proxy is not in the right place, but I don't think the first problem (without problem) comes from here.


texts: Displaying quantitative data in choropleth legend?

texts:
I have always struggled to come up with a great solution to handle choropleth or any other quantitative interval or dataset classifications in a map legend where the numbers are sequential.

For instance: here is the ArcGIS default legend, which is a bit gaudy and what's the real cuttoff? .55000000001 or .5501 or what?:

And I see many people create their legends like this:

Usually with "TIE" written as "50%", which in this case would be wrong on two levels, 1. These election results (like many) have pluralitiy winners,i.e. winning with 49.8% 2. If a county/unit is 55%, which classification is the unit a member of? It can't be both! And this particular dataset does contain 4 tied units.

I tend to use interval notation in my legends, like this:

Set notation { x| 55 < x <= 67} would work as well, but I'm worried that people won't always understand these methods.

Anyway, I'd like to hear your solutions to this to help figure this out.


texts: QGIS 3.18 Override raster gradient legends

texts:
In the new QGIS 3.18 (Zurich), raster color ramps automatically show up in the print layout as gradient legends! This is great, but I am trying to match some old maps to new maps for a publication and would prefer the old categorical legend type.
The change log mentions there is a way to override this behavior, however I have not been able to find it.
(https://www.qgis.org/en/site/forusers/visualchangelog318/index.html#id32)
Any ideas?


texts: What is the order of points of a polygonz surface in postgis

texts:
I extracted points of polygonz surfaces using ST_DumpPoints. I need to generate grid points on these surfaces so, for this purpose I need to know order of points. Because I need to know which vertices are neighbour and which are not. I am using this formula to calculate coordinates of new points ;

Every piece has base vertex A0 (x0, y0, z0) and two neighbor vertices A1 and A2. Make base vectors

V1 = A1 - A0 V2 = A2 - A0 To make regular grid points inside this piece of plane,

P(u, v) = A0 + V1 * u + V2 * v where u and v are equally-spaced parameters

So what are the order of points that make polygonz surface?


texts: Select line on the map and check its attributes in QGIS

texts:
In QGIS can I select a line in the map and check its attributes?

I have an ESRI shape file but I am not able to select or choose any line of the map in QGIS, I see, it looks like one block but in show attribute table, I see lots of different id.

When I show attributes in the table and select any id, QGIS should show me which line belongs to this id on the map, right?

I am totally new and I need to have good start in it.


texts: get only property name and value from getGetFeatureInfoUrl WMS layer ol5

texts:
Stuck in a problem to get only property name and value of WMS layer.
here is my code

var point_wms = new ol.layer.Image({
  title: 'Gound Control Point',

    source: new ol.source.ImageWMS({

        url: 'http://localhost:8080/geoserver/wms',

        params: {
            'LAYERS': 'export:Ground_Control_Point',

            'VERSION': '1.1.1',
            'FORMAT': 'image/png',
                'tiled': true,
                serverType: 'geoserver',
                crossOrigin: 'anonymous'
        }
    }),
     visible: true,
    opacity: 1
});

map.on('singleclick', getFeatureInformation);

function getFeatureInformation(evt){
        var coordinate = evt.coordinate;
        var lonLat = ol.proj.transform(coordinate,'EPSG:3857','EPSG:4326');
        var viewResolution = (map.getView().getResolution());

        if(point_wms.getVisible()){
            var url = point_wms.getSource().getGetFeatureInfoUrl(
              evt.coordinate, viewResolution, 'EPSG:3857',
              {'INFO_FORMAT': 'text/html'});
            console.log("url = "+url);
                if (url) {
                  content.innerHTML ='';
                }

                 overlay.setPosition(coordinate);

            }
}

it shows result like

but i want to get property name and value so that i can show this result using my html code in desired format like vertical table. how to achieve this?

i edited the question here is my code

 if(point_wms.getVisible()){
            var url = point_wms.getSource().getGetFeatureInfoUrl(
              evt.coordinate, viewResolution, 'EPSG:3857',
              {'INFO_FORMAT': 'application/json'});
                /*console.log("url = "+url);
                if (url) {
                  content.innerHTML ='';
                }*/
            if (url) {
                var parser = new ol.format.GeoJSON();
                console.log(parser)
                $.ajax({
                      url: url,
                      dataType: 'jsonp',
                      jsonpCallback: 'parseResponse'
                    }).then(function(response) {
                      var result = parser.readFeatures(response);
                       console.log("result.length = "+result.length);
                      if (result.length) {

                        var iso3 = [];
                        for (var i = 0, ii = result.length; i < ii; ++i) {
                          iso3.push(result[i].get('Point_Id'));
                        }
                         content.innerHTML = iso3.join(',');
                         overlay.setPosition(coordinate);
                      }
                    });
              }

            }

still not getting value of point_id.

it is showing error like this

Uncaught SyntaxError: Unexpected token <
    at m (jquery.min.js:2)
    at Function.globalEval (jquery.min.js:2)
    at text script (jquery.min.js:2)
    at Ut (jquery.min.js:2)
    at k (jquery.min.js:2)
    at XMLHttpRequest.<anonymous> (jquery.min.js:2)

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://localhost:8080/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException code="ForbiddenFormat">
      Getting feature info using application/json is not allowed
</ServiceException></ServiceExceptionReport>


texts: What is the difference between GRASS 6.4.0 and GRASS 6.4.SVN?


texts: qgis2threejs export quality of map canvas image

texts:
A comparison

qgis2threejs preview: quality near the original within qgis: Very good

browsers chrome+edge: bad (edge a little bit better than chrome)

firefox: extremly bad

What can I do to improve the quality within the browser?
Here the four screenshots as one image:


texts:
Can you have an ArcGIS Online map, you know one that you share and embed on a web page or link to, automatically start at an extent that is defined in a query of some sort?

Ideally it would be great if it could be inserted into the HTML link such as:

https://domain.maps.arcgis.com/apps/webappviewer/index.html?id=a618e74cc2cb439197413549eeb1256bx123&Fieldname=XYZ

Can this be accomplished or something similar?


texts: Applying color style to GeoTIFF layer using Mapnik and Python

texts:
I am trying to apply a color style to a GeoTIFF using Mapnik.
The code I'm using is:
import mapnik

map = mapnik.Map(800,600)
map.srs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'

style = mapnik.Style()
rule = mapnik.Rule()
sym = mapnik.RasterSymbolizer()

rc = mapnik.RasterColorizer()
rc.colorizer = mapnik.RasterColorizer(mapnik.COLORIZER_LINEAR, mapnik.Color("white"))
rc.colorizer.add_stop(-20, mapnik.Color("gray"))
rc.colorizer.add_stop(0, mapnik.Color("yellow"))
rc.colorizer.add_stop(10, mapnik.Color("red"))
rc.colorizer.add_stop(15, mapnik.Color("green"))
rc.colorizer.add_stop(25, mapnik.Color("blue"))

sym.colorizer = rc
rule.symbols.append(sym)
style.rules.append(rule)

map.append_style('style', style)

layer = mapnik.Layer('My data', map.srs)
layer.datasource = mapnik.Gdal(file='input.tif')

layer.styles.append('style')

map.layers.append(layer)

mapnik.render_to_file(map, 'map-output.png', 'png')

The GeoTIFF gdalinfo output is:
$gdalinfo input.tif
Driver: GTiff/GeoTIFF
Files: input.tif
Size is 1414, 631
Coordinate System 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","4326"]]
Origin = (-125.229333333000000,50.146000000000001)
Pixel Size = (0.041700000000000,-0.041700000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-125.2293333,  50.1460000) (125d13'45.60"W, 50d 8'45.60"N)
Lower Left  (-125.2293333,  23.8333000) (125d13'45.60"W, 23d49'59.88"N)
Upper Right ( -66.2655333,  50.1460000) ( 66d15'55.92"W, 50d 8'45.60"N)
Lower Right ( -66.2655333,  23.8333000) ( 66d15'55.92"W, 23d49'59.88"N)
Center      ( -95.7474333,  36.9896500) ( 95d44'50.76"W, 36d59'22.74"N)
Band 1 Block=1414x1 Type=Float32, ColorInterp=Gray
  Min=-23.065 Max=24.441
  Minimum=-23.065, Maximum=24.441, Mean=0.561, StdDev=6.495
  NoData Value=-3.39999999999999996e+38
  Metadata:
    STATISTICS_MAXIMUM=24.44076094814
    STATISTICS_MEAN=0.5611245629664
    STATISTICS_MINIMUM=-23.065323150453
    STATISTICS_STDDEV=6.4950334131849

The problem is that the PNG result does not show the colors I'm applying. It shows only a empty PNG image.


texts: Process finished with exit code -1073741819 (0xC0000005) from RegionGroup (ArcGIS 10.7)

texts:
I'm trying to automate my work using ArcPy in PyCharm on Windows. However, when I run the following (abridged) loop, it exits at random points: sometimes at the 4th iteration, sometimes at the 31st, sometimes it manages to get to 100 but never close to the 90.000 iterations it's supposed to!

The error is:

  Process finished with exit code -1073741819 (0xC0000005)

import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
input_path =  "D:/.../05_reprojected"
current_path = "C:/.../Arcpy-test"

for filename in os.listdir(input_path):
    glacierID = glacierLocation[27:35]+"_"+glacierLocation[36:41]
    bedrockRaster = Raster(glacierLocation)
    filledBedrockRaster = Fill(bedrockRaster)
    fill_save = current_path+"/filled_"+glacierID+".tif"
    filledBedrockRaster.save(fill_save)
    sinks = Con(fill_save > bedrockRaster, 1)
    regionGrp = RegionGroup(sinks)
    # sometimes it breaks here, and sometimes it doesn't
    region_save = current_path+"/regionGrp_"+glacierID+".tif"
    regionGrp.save(region_save)
    ...
    # and the loop continues

I know from googling that this has something to do with Windows access violation errors. The regionGroup(sinks)output in the terminal is:

Starting region group
Connecting regions
Compacting labels
# and here it breaks sometimes with the above mentioned error

Where would this access violation take place?

I store every output in a separate folder for every iteration. If RegionGroupcreates an own tmp-folder, I don't know how to find it and check for issues.

I have looked at many questions on GIS Stack Exchange. There is an ArcGIS (Desktop, Engine) Background Geoprocessing 64-bit Arcpy Exit and Shutdown Patch mentioned at ArcPy 64bit giving 0xC0000005 error with PyCharm, but it is not for ArcGIS Desktop 10.7, which seems to indicate that it should be fixed in this version (but is not, apparently). Other questions are related to PyQGIS or other workflows that do not apply to mine.

EDIT:
I've tried to run the script via the ArcGIS python directly in the console and but the error was the same. So it does not seem to be connected to PyCharm.


texts: Getting reference to existing reference of dockable window in class in ArcMap Addin?

texts:
I have an addin for ArcMap 10.2 that has a dockable window.  For what it's worth, this dockable window has a DataGridView to which the user can add features for review.  These features are added from other buttons and menus that are part of the addin.  All of that works fine, as long as the dockable window is not showing when the MXD is loaded.

However, if the dockable window is showing at startup (for example, if the user left it open and saved their MXD at the end of the last session), then I am unable to access the dockable window (or add data to the dockable window), and I get a Null Reference Exception.  I assume that there is an already existing copy of the dockable window class somewhere (since the dockable window is visible in ArcMap), but I'm at a loss on how to access it.

I've tried using the DockableWindowManager to get a reference to the dockable window and then close it; I've also tried to see if an existing copy of the class can be found using Marshal.GetActiveObject.  Neither of these solve my problem.

Is there a way to get a reference to an existing copy of a dockable window class that I can use?


texts: Define your own latitude and longitude projection

texts:
I am converting (rotating) vectors from Lambert to latlon (and viceversa) in Python. But in Basemap or Cartopy, I see that you do not define the lat and lon projections, you just define your Lambert, but the system always take lat and lon as (I guess) the WGS-84 datum. For example, from this example here (https://stackoverflow.com/questions/8878564/how-can-i-rotate-vectors-onto-a-basemap-map-projection):
import numpy
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

# Set up map projection
m = Basemap(llcrnrlon=-15.,llcrnrlat=46.,urcrnrlon=15.,urcrnrlat=59.,
        projection='lcc',lat_1=40.,lat_2=50.,lon_0=-50.,
        resolution ='l')

# Calculate positions of vectors on map projection
x,y = m(lon,lat)

How can I define my own latitude, longitude system (I want a spherical one, not the ellipsoid of WGS-84) in order to perform the conversions from my Lambert to my latlon? I can use Cartopy, Proj4 or whatever Python tool you suggest.


texts: Arcobjects 10.3 C#: How to re activate the datasource of an ILayer, in memory

texts:
I captured (stopped) a Layer (featureclass) from being deleted from the Table of contents of arcmap. When i re-add it back, the datasource is broken.

(See question :Using ArcObjects to stop items from being deleted from Table of Contents?)

I'm wondering what I need to cast my ILayer object to, in order to just activate the datasource?

The workspace is sde, but i find it odd that if the whole object(layer) is available when caught by the itemdeleted handlers, why is the datasource off?

I'm having trouble finding the right interface to use.

Here is the code that I have implemented: It readds the layer back to the ToC but the connection is not connected.
(Exten is the Addin.extension object for the addin.)

 public void On_TOCItemDeleted(object Item)
        {
            IWorkspaceName workspace;
            var msgText = "You tried to delete an item from the TOC.";

        // check that the item is of type ILayer (lyr file, raster, etc) before doing anything.
        var DeletedLayer = Item as ILayer;

        if (DeletedLayer == null)
            throw new Exception("Not iLayer file");

        IDataLayer2 src = DeletedLayer as IDataLayer2;
        var groupname = FindParentGroup(DeletedLayer.Name);

        if (src.DataSourceName != null)
        {
            workspace = (src.DataSourceName as IDatasetName).WorkspaceName;
        }
        else
        {
            switch (groupname)
            {
                case BLayerName: // basemaps etc
                    workspace = DeletedLayer.SetLayerSDEConnection(Exten.Analysis_WIP);
                    break;

                default: //analysis connection
                    workspace = DeletedLayer.SetAnalysisSDEConnection(Exten.Analysis_WIP);
                    break;
            }
        }
        src.Connect((workspace as IName));
        msgText = string.Format("You tried to delete the layer '{0}'. This operation was stopped.", DeletedLayer.Name);
        //Read the item

        var grouplyr = GetGroupByName(ArcMap.Document.FocusMap.Layers, groupname);

        grouplyr.Add(src as ILayer);

        ArcMap.Document.UpdateContents();
        ArcMap.Document.ActiveView.Refresh();

        Exten.InformUser("App", msgText);
    }


texts:
How I can get Digital Height Model (only object heights, DSM - DTM = DHM)?

Many applications can do this, but they convert it to GRID or TIN format. I was wondering if it was possible to keep the data as point clouds.


texts: QGIS / QField: nested relational value form QGIS save refused QField

texts:
The Field Survey project should have been operational for two days already, but I have a major problem.
After having succeeded in nested two fields in the form (relational value), we realize that the saving under QField is not done!
Where is the concern?
QGIS side, in the form? (group/detail interlocking is not active in update mode)
QField side? who if he does not know how to read the expression, does not want to take it into account?
I have attached several screenshots to understand my concern
Hope any of you can help me solve this blocking problem.
Additional info: I switch from QQGIS to QField by QField sync.


texts: I am looking for a free WMS layer with road map. It needs to support NAD27 so we can display some old data on top of it. Long time ago we used to work with google maps; but they suspended the service unfortunately. I tried OpenStreetMap; but I did not see support for NAD27 (if I'm wrong let me know).


texts: Accessing ArcGIS Online feature dataset from R using R-ArcGIS Bridge

texts:
I can successfully connect to an ArcGIS Pro feature dataset stored on my local computer in a file geodatabase in R like so:

library(arcgisbinding)
library(sf)
arc.check_product()

districts <- arc.data2sf(arc.select(arc.open("Maps/MyProject.gdb/districts")))

I'm trying to figure out whether it is possible to connect to files that are stored and updated on ArcGIS Online. These files are not "public", they are only available to myself and my team. Essentially, features that are part of "My Content" online.

I've tried the following using the URL available on the feature's page, but no luck (this is not a 'real' URL):

test <- arc.data2sf(arc.select(arc.open(
"https://services1.arcgis.com/jsdf8234sdfs/arcgis/rest/services/districts/FeatureServer")))

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘arc.select’ for signature ‘"arc.workspace_impl"’


texts: Remove small areas from SHP shapefile in Qgis? (Python R etc)

texts:
I have a very complex file showing remnant vegetation data. I wish to remove very small areas (preferably with a threshold that I can set) and generate a file that shows only larger blocks of vegetation.

Is this possible with Qgis (2.14.3 Essen) on Ubuntu. I could also do this in Python or R if someone can give me a clear script (I can just about hack R).


texts:
I am trying to render the MapQuest Open style at zoom level 3 for the whole world using Mapnik, at a higher resolution than available online.

I have downloaded mapnik, and the MapQuest Open style. However, downloading the planet.osm and importing data into postgis is beyond my capacity and seems unnecessary for such a low zoom level. All I really need is place names, borders and land cover data. I have downloaded some shapefiles containing some of this data, however I do not know how I would import this into postgis or how to use these with the MapQuest Mapnik Style.

Is it possible to render a very low zoom level map of the world with mapnik without downloading large amounts of data? If so what is the easiest way to go about this?


texts: Python Program for "QC of PDFs Created from Data Driven Pages"

texts:
I need some help. I am fairly new to Python. I found a sample code and I need to change a couple things in it for it to work. I need it to make sure that all PDFs were created from the Data Driven Pages program I have now. Reason is I sometimes have about 150 PDFs created from the DDP program. For Example : Houston has 152 pages in my mxd and DDP is supposed to make 152 pages. I've made some changes already to the code I was needing someone to double check my work.

import arcpy, os, sys
from os import path as p

#Make parameters for people to choose mxd and folder of PDFs to compare
mxdList = arcpy.GetParameterAsText(0).split(';')
pdf_path = arcpy.GetParameterAsText(1)
fc_path = arcpy.GetParameterAsText(2)

report = {}
page_list = []

#Use Search Cursor to go through Attribuite Table to get Sheet number info
for mapDoc in mxdList:
    mxd = arcpy.mapping.MapDocument(mapDoc)
for fcClass in fc_path:
    fc = arcpy.ListFeatureClasses()
            max_list = []
            rows = arcpy.SearchCursor(mxdList, "Sheet_ID")
            for row in rows:
                max_list.append(row.Sheet_ID)
                sheet_count = max(max_list)
                arcpy.AddMessage('Sheet count: {0}'.format(sheet_count))
            else:
                arcpy.AddError('No Layers in %s match data source'%mapDoc)
pdfDoc = arcpy.mapping.PDFDocumentOpen(path)
print pdfDoc.pageCount


texts: What is the difference between ArcGIS Desktop and ArcGIS FOR Desktop?

texts: I keep seeing these two categories of software package on the ESRI site but nothing making it clear the difference between what ArcGIS Desktop is vs. ArcGIS For Desktop.  I use ArcGIS Desktop (i think) that has ArcMap, ArcCatalog, etc. what does ArcGIS For Desktop have that is different?


texts: How to save .shp with layers 'baked' in with QGIS

texts:
I'm having no problem applying labels to the shapefile I'm working on in QGIS, but I can't seem to save the new shapefile with the layers 'baked in' for uploading to Mapbox, and then styling.

I've tried lots of saving options after I apply the labels, but can't seem to get them to show up when I import it into Mapbox Studio. Any help?


texts: Check multiple feature classes have same spatial reference

texts: As the title says I've several feature classes. I want to check whether they have same spatial reference using arcpy


texts: selectFeatures Javascript Removes Feature

texts:
I've been messing around with linking up a dgrid to select a feature layer and then fire an info template.

When it fires, it removes the selected point feature.  Right now I have it set to when you click on the objectID in the 'ID' field of the dgrid it will select the feature.

Can view here on fiddle: http://jsfiddle.net/t2xan9pL/


texts: Making subsections in network-dataset for using location allocation in ArcMap?

texts: I'm making a simple smart grid model on ArcGIS for my thesis and I have a town surrounded by several renewable energy sources from all around. I want to divide the network dataset into four subsections so certain energy sources supply certain subsections (which will each have a demand point). I'm using ArcGIS (ArcMap).


texts: Creating graduated point circles in ArcGIS Desktop?

texts:
I have not been giving any data, but I have been told that I could be expected to recreate maps In ArcGIS 10. Below is an sample portion of what I may have to recreate. I have no idea  how to get symbology like this for point data. Graduated symbology would only work across the array of points. How do I get one point to have the color variations like seen below? The points are representing observation density. I know how to obtain and symbolize basic data. I don't need information on how to get the hydrology layers or symbolize the eagle point. I need to know how to do the red circles.

Heat maps are not what I need since the "blanket" the entire area.


texts: Combining rasters after weighted overlay

texts: I have performed five weighted overlay analyses in arcgis 10.2 and got five rasters with pixel values ranging from 1 to 10 as outputs. Now I need to combine those rasters in one single raster stating layer name which is having maximum dominance for each pixel. Cell statistics will only give the actual maximum pixel value among several input rasters, but I need to identify raster 'file name' which is having maximum value, for each pixel.


texts:
I'm trying to calculate geometry in Python for a script tool for ArcMap.

For the script I want to calculate area in hectares. I am trying to use the add geometry attribute (data management) tool, but it's not working and it's giving me all zeros in my output field. This is what I have:

# Set Calculate geometry variables

calculate_table = clip_output

geometry_prop = "AREA"

expression = "HECTARES"

# Run calculate geometry

arcpy.AddGeometryAttributes_management(calculate_table, geometry_prop, "", expression, "")

What's wrong?
Should I use a different tool?


texts: Serving features of mixed transparency levels from ArcGIS Server map service?

texts:
I'm using a dynamic map service layer that is shown in a semi transparent way overlaid on a tiled map service.

There are a few polygon layers in the dynamic map service layer that cover the whole mapped area in different solid colors. Those layers are set as transparent in the MXD, but when the map is served out from ArcGIS server, it displays differently from in ArcMap. The layers blend down with the background MXD itself, making them turn %50 white, and making the basemap behind them as difficult to read as if they were solid.

Is there a way to serve features of mixed transparency levels from an arcgis server service?


texts: Interactive "for each item in ... do ..." in QGIS?

texts:
How does one interactively run a basic looping batch construct in QGIS?

for each {feature} in {features}:
    do {thing} with/to {feature}
    save/export/display {result}

Where:

{features} could be a polygon shapefile, a table, a selected record set, ...
{feature} an individual record from {features}
{thing} is any Processing tool (or menu action?)

In ArcGIS cursors, batch processing, and model iterators can be brought to bear for looping, though only the latter are interactive. What are the QGIS equivalents?

Most of the topics and tutorials I've looked at so far involve shell scripts or straight up python code -- meaning little to no interactivity and dynamism. Is this really the best/only way, at present? Essentially, I'm re-asking Interactive "For each .... Do" loop? and hoping for a QGIS specific answer.


texts: GDALwarp Creates Darker target Image sometimes

texts:
I'm yoyo and new to this forum. I'm relatively new to GIS and GDAL. Please accept some typos now and then... English is not my native language.

Because of some problems with GDAL 1.7.0b2 in FWTools 2.4.7 I upgraded to GDAL 1.11.0 in the latest OSGeo4W distro. This solved some issues I had, but introduced a really anoying one and what ever I do I can't get rid off it. I hope someone can help me out. If btw I try to go back to GDAL 1.7.0b2 with FWTools the problem is still there, so I guess it's a driver/DLL update issue.

My problem looks exactly like this post... GDAL TIF to JPG Creates Dark Image
Only that I don't use JPG. And I use GDALwarp. And not all target files have the problem. Most target files look ok.
And I only see it with Aerial Photo like sources. With vector like raster sources I don't see the problem. So when reading the other post I thought it had to do with a mask/alpha channel or something, but I can't see this in the source. And whatever option I test, the problem remains.

Section of Source : https://www.dropbox.com/s/gy2tu0h8dufv41o/Source.jpg?dl=0

Intersection of four Target files : https://www.dropbox.com/s/av7mdj8ivx3zgsl/Target.jpg?dl=0

I tested almost all options in multiple actions in my workflow.
So cutting the source with gdaltranslate in smaller parts for editing. Building a VRT with gdalbuildvrt and cutting into final parts with gdalwarp which will show the problem.

I started with the -b 1 -b 2 -b 3 option in my gdalbuildvrt action, but GDAL crashes here... (?)
I also tested other options that might have something to do with the problem...
-nomd -co ALPHA=YES/No -co PHOTOMETRIC=RGBA/RGB -setci, you name it...

I checked good and bad source and target areas with gdalcompare and gdalinfo, but I saw nothing strange...

I'm now thinking of resetting GDAL driver settings if that is possible..(?)

Has anyone seen this before? Can anyone give me a hint..?

PROBLEM FOUND.

I'm using a workaround now...
I found out that the problem was there with one specific target size. So I started testing with filters and when I used an other filter then cubicspline, the problem was gone. So the problem was a combination with the cubicspline filter and a specific target resolution.


texts: Determine distance between two polygon with having unique id

texts: I have two polygon feature class one is Source Data and other is shifted data. I want to calculate how much polygon shifted (distance) using ArcGIS 10.6.


texts: How to import USGS SDTS ddf file into Qgis

texts:
I have downloaded SDTS dem files from the USGS.  These extract from tar.gz files as ddf files.
I have tried importing these files into Qgis using "add raster layer and choosing .ddf or .DDF from the file type list.  But the program returns "Unsupported Data Source: C:\tmp\RockyHills_dem\1694744.DEM.SDTS\9948CEL0.DDF is not a supported raster data source".

Can anyone help?


texts: Using GPU for processing Inverse Distance Weighted within Project Modeler

texts:
We are currently doing IDW interpolation within QGIS’s Project Modeler only it is taking over 24 hours to do one run.  And we have over 30 to do.
Is there any way we can get it to process the data using the computer’s GPU’s rather than the CPU in order to (significantly) increase the processing time?
I have heard about using Python’s "Multi Processing" however do not know how this could be integrated into the Project Modeler.


texts:
I am working in a tabular set (.dbf) In ArcGIS 10.1 with an odd DMS format.The format is as follows:

Lat: 45585400
Long: 89560600

I believe this is a strange notation for 45°58'54" and -89°56'06". My problem is I can't figure out how to parse the Degrees, Minutes and Seconds out into individual fields in the .dbf. Once I get that far I can easily convert the value to decimal degrees and create point data out of it.

If anyone can help me with this I would really appreciate it.


texts: Newbie help -- What do I need to calculate the 3D area of a terrain?

texts:
I have a terrain right on the side of a hill in Mexico. Initially I thought I could get the area from Google Earth, but after a few tests I found out that the area Google Earth calculates is flat, i.e. it does not take into account the undulation and topography of the terrain. This is a problem, because for legal purposes I must find out the actual area of my terrain taking into account the terrain's topography.

Now, I do not know anything about GIS, but digging around Google turned up what looks like a surface volume module of 3D Analyst for ArcGIS. I intend to locate a civil engineer or an expert in geography and have them do it. How do I ask them to do it, and what information if required must I provide to them in order to do that?

The module in question is this one: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00q900000027000000


texts: WHERE clause that's stored in a table: Use it to filter a layer?

texts:
I have a workorder management system (WMS).

The WMS has functionality that allows users to save queries.
The queries are saved in a Oracle database table called QUERY (in a CLAUSE column).

An example of a value in the CLAUSE column would be:

((woclass = 'WORKORDER' or woclass = 'ACTIVITY') and contains(description,'$paint') > 0  and historyflag = 0 and istask = 0 and siteid = 'SERVICES')

I would like to utilize the saved WMS queries by applying them as filters to layers in a Portal web map.

From a drop-down menu, a user would choose a CLAUSENAME (from the QUERY table), the widget would programmatically extract the WHERE clause value from the CLAUSE column, and use it to automatically filter a layer.

To be clear, this would be a programmatic/custom widget solution, not manual copying and pasting.

Is it possible to do this?


texts: Measuring river width every mile in Google Earth?

texts:
Given an area, can Google Earth  tabulate the width of a particular river?

This problem I think requires using an API to programmatically:

define an area of interest
identify a terrain feature like a river
iterate over the length of that terrain feature measuring width

Although it sounds simple and a child could do it using the google earth application and the ruler tool, it seems hard to write code that could do this.

A simpler problem (still hard) is measuring the length of a river and its tributaries in a defined area. Is there a tool in Google earth that can do this?

I realize their API has been deprecated so I am open to another application/solution for this.


texts: Working around Table does not contain OIDs arcpy Error 000339?

texts:
I have a script I'm working on that takes a SQL table view and performs a Table to Table (arcpy.TableToTable_conversion)

If I leave out the Where Clause the Table to Table works flawlessly, however if I include the Where Clause it fails throwing Error 000339 "Input does not have OIDs".  Typically the way to overcome this error is to copy the table (or table to table) but this is what I'm doing.  The table is quite large and takes some time to copy, which is why I'm trying to add the Where Clause.

Is there another way I can get my table copied, filtered by Where Clause, that won't throw the missing OID error?

whereclause = "FeatID IN ('501684', '179845', '173458', '509634')"
arcpy.TableToTable_conversion(WMTR_Pivot, TargetDB, "xxTest_WMTR_Pivot", whereclause, wmtr_fms, "")

Note: FeatID is a text field.

Error:

Traceback (most recent call last):   File "xxAppendMeters.py", line 231, in <module>
    arcpy.TableToTable_conversion(WMTR_Pivot, TargetDB, "xxTest_WMTR_Pivot", whereclause, wmtr_fms, "")   File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\conversion.py", line 2133, in TableToTable
    raise e ExecuteError: ERROR 000339: Input v_WMTR_Pivot does not have OIDs Failed to execute (TableToTable).

My table view has about 75000 records, and I'm trying to filter that to about 150 records with the Where Clause.  I've trimmed that to 4 in the code snippet above.

Table is in SQL Server 2012 Enterprise Geodatabase, on my PC I'm running ArcGIS Desktop 10.3.1


texts: Filter OGR datasource in Python (or XML)

texts:
I'm trying to streamline the process whereby we build tiles (of ~30 layers, representing planning zone and zone-overlay boundaries). I have a strong preference for FOSS tools.

A couple of our people - and a couple of key clients - are only really comfortable in MapInfo, but I'm comfortable converting across platforms (using OGR) because I prefer PostGIS/QGIS/Python. However the source files are always (or nearly always) MapInfo TAB.

Many of the TAB files are 'thematic'; the theme style results either from the values taken by one aspatial field (e.g., 'LGA'), or from the interaction of 2 aspatial fields (e.g., 'Zone' and 'Status'). The thematically-different combinations ought to be different layers (or better yet, different files) but they ain't: I have no control over that, otherwise it would change and this question would be moot.

So I'm trying to separate out these layers, style them (that's easy) and then chop them into a maptile pyramid (super-easy). The eventual aim is to build XML files so that (python + mapnik) can style each layer, then chew it into a tile pyramid.

I'ma do it in Python with mapnik and ogr.

Now

ds = mapnik.Ogr(base=[dir], file =[tabfile], layer_by_index=0)`

or just

ds = ogr.Open('[tabfile]')

brings a TAB file into a Python script without any problem - analogous entries in the <Datasource> tag in an XML file would do likewise.

In the mapnik version, ds.fields() shows me the right headers, and I can add rule filters to extract specific subsets of the data to style as I wish, e.g.,

m=mapnik.Map(800,600)
s=mapnik.Style()
r=mapnik.Rule()
r.filter = mapnik.Expression("([Zone]='IN2Z') and ([Status]='n')")
s.rules.append(r)
m.append_style('IND_01', s)

What I would like to do, is to loop through distinct values of, say, Zone and Status, and filter a new dataset (and apply a new style) at each step.

If there was a way to return a vector of the unique values of the relevant field, I could simply iterate through the vector and do my dastardly work: for example, if I have a TAB file BH.TAB with column 'LGA', at the command-line I can do

ogrinfo -sql "select distinct LGA from BH" BH.TAB

This (piped to a file, say) can be parsed to fetch the vector of unique LGA in BH.TAB.

So now - finally - the question.

Does anyone know ogr's Python bindings well enough to know how to emulate the ogrinfo command-line syntax shown above? Once I can get that sucker into the mix I will find a way to parse it.

I swear I have RTFD (e.g., the Python/OGR Cookbook, which is pretty comprehensive), to no avail.

UPDATE:

I should point out that I have already tried stuff like

v = subprocess.check_output("ogrinfo -sql \"select distinct LGA from BH\" BH.TAB")
w = v.split("\r\n")

The output of that is a PITA to parse.


texts: Changing data source on all layer files in folder using ArcPy?

texts:
I have a folder of .lyr files and I want to change the source on all of them from one file geodatabase to a different file geodatabase.  Feature dataset and feature class names will remain the same in both geodatabases.  There are a few layerfiles that are grouped layers and there may be times when the FC doesn't exist when the script is run but will be there in the future which is why I have validated False.  I'm a beginner at using Python but this is what I've managed to throw together based on my research.

Using the Python window in ArcMap I've managed to change data source on the .lyr files that are not group layers with the code below, but I cannot get it to work in a .pyt toolbox where I hoped to set up parameters for the Layer File folder and the new geodatabase.  The script runs in the .pyt tool box but I get no messages and the layerfiles are not updated.

I also cannot figure out what I need to do to change the data source on the group layers.

Not sure if it matters but this is on ArcMap 10.5.1 run in a Citrix environment.

For the grouped layers:
I've tried using glyr = arcpy.mapping.Layer(lyr) with error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 8, in <module>
  File "e:\arcgis\desktop10.5\arcpy\arcpy\arcobjects\mixins.py", line 428, in __init__
    raise TypeError(repr(type(lyrfile)))
TypeError: <class 'arcpy._mapping.Layer'>

and glyr = arcpy.mapping.ListLayers(lyr)  with error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 9, in <module>
AttributeError: 'list' object has no attribute 'name'

Python code that works in the ArcMap Python window:

import os
for dirname, dirnames, filenames in os.walk(r"T:\Project\TestLayerFiles"):
    for filename in filenames:
        newFileName = os.path.join(dirname, filename)
        if not newFileName.lower().endswith(".lyr"):
            continue
        lyr = arcpy.mapping.Layer(newFileName)
        if lyr.isGroupLayer:
            continue
        if lyr.supports("WORKSPACEPATH") == True:
            print("DS :" + lyr.dataSource)
                lyr.findAndReplaceWorkspacePath(r"T:\North\ProjectA\projectA.gdb", r"T:\North\ProjectB\projectB.gdb",{False})


texts: OpenLayers 4 and MapServer 7.0.4 HTML template

texts:
I would like to know if it is possible to render directly a HTML template retrieved from a MapServer WMS GetFeatureInfo call.

The documentation explains that the API uses ol.Format.GML2 to read features Here the API DOC. No mention about HTML.

Here the template I would like to render: HTML template by MapServer

I would like to not know if there is a "direct" way, not as shown in this
example.


texts: Emergency Dept. Distance OSM roads Network Analyst

texts:
I have a shapefile that is all of the roads in the Republic of Ireland, I have downloaded this from OSM. I have then created a new point shapefile with a point for every emergency department in Ireland. I am trying to find out what road or section of road is the furthest from an emergency dept. in the country.

From what I know, I can clean up the geometry and then create a network dataset and then maybe use 'Closest Facility' to find what I'm looking for?

Will I be able to create a network dataset from ordinary OSM data?


texts: "Select layer by location" tool - input feature layer Error

texts:
I am trying to use  the "Select layer by location" tool, in Model Builder, ArcGis.
This tools has 2 input files: "Input Feature Layer" and "Selecting Features".
I cannot connect the "Input Feature Layer".
When I drag and drop the "Connect" arrow (which is used to link layers to tools), The tool only shows "Selecting Features" (Fig. 1). Ok... I have tried to open the tool, and choose from there, the "Input Feature Layer", but it creates a duplicate of that layer (Fig. 2). By creating a duplicate, all the steps of this model stop at this particular step.
What can be causing this error? I work in ArcGis 10.2.2 and use no extensions.


texts: Migrating Leaderline Annotations from AutoCAD to GDB using FME?

texts:
I am trying to migrate leaderline annotation from AutoCAD to ESRI thru FME. While reading the leader anno from dwg, the reader reads it as 3 individual components and not a single entity.

Is it necessary to define all geodb_text_leader* attributes?

Currently, only the text gets migrated and not the leader line.


texts: Geoprocessing tool for finding specific point within polygon

texts:
I have a ArcGIS task and cannot find a solution.
I try to find a certain point within a polygon (see picture). This example is handmade and I need a geoprocessing solution

I point I try to get is the one that is furthest away from the borders of the polygon, but not only from the outline of the polygon (this would be centroid), but ALSO from the inner line. Thus this would be the point around which I could draw/buffer the biggest circle around within the polygon (without touching the borders of the polygon).
I hope my explanation was understandable.


texts: Administering 10.3.1 ArcGIS Server services using ArcGIS Server 10.3.1 and 10.2.1

texts:
I am using the following line of code to stop ArcGIS Server services in a 10.3.1 version.

C:\Python27\ArcGISx6410.3\python.exe "C:\Program Files\ArcGIS\Server\tools\admin\managesite.py" -u <USERNAME> -p <PASSWORD> -s https://<SERVERNAME>:6443 -c default -o stop

The code works fine when run on the ArcGIS Server 10.3.1 machine.

Is it possible to stop these 10.3.1 services (using this line of code) from another machine where ArcGIS 10.2.1 is installed? That is, can I run this code from the 10.2.1 ArcGIS Server machine to stop the 10.3.1 services?


texts: Converting Spidertracks GPS flight tracks into shapefile or kml?

texts: I have some data that will be comming to me from helicopter flights and was recorded with Spidertracks software. Does anyone know how to convert flight tracks from 'Spidertracks' gps flight tracks into shapefile or kml? What is the raw spidertracks gps format? I am trying to learn before I recieve the data and I'm not 100% sure what format I will be recieving.


texts:
How I could extract a cube root of the multiplication of 3 rasters in raster calculator?

Eg: cuberoot (ras1 * ras2 * ras3)


texts: How to keep NA values when using st_intersection in R

texts: I have a polygon layer representing different site treatments and a point layer representing plots taken within those treatments and several control plots outside of the treatments. I need to extract the treatment information from the polygons and add it to the point layer, but there are several control plots that are outside of the polygons that I need to keep. I have been using the sf package in R. Is there a way to use st_intersection and keep NA values for the points that don't intersect the polygon layer, or do I need a different strategy?


texts: HEATMAP-Nan Layer, Kernel, Rows and Column

texts:
My Friends,

My problem is with the Heatmap Plugin.

I have a 20.000 points in the map and when i want to heatmaping, the layer is NaN.

the village has 88 Km2 and is in the OpenStreetMap (wgs84 epsg:4326)

The values of heatmap are:

Radius: 300 meters

Rows: 50000
Column: 77712
*i prove with 500, 5000 and 50000 and i have the same problem (nan)
xcell: 0.00187496
ycell: 0.00187496

If someone can help me, please!


texts: Using US Census block group data field names for feature classes?

texts:
I may be missing something here, but it seems like getting usable US Census data into a shapefile or feature class is much harder than it used to be. Downloading the Tiger block group shapefiles is no problem, but getting and joining the data is.
It used to be you could just drill down to the table and area you wanted, download, alter an ID field so you could join data, and you'd be done. Now it seems hideously complicated.
As a result, I downloaded the Tiger gdb that the census says is

  A limited set of TIGER/Line Shapefiles are available pre-joined with data in geodatabase and shapefile format.

Ok, block groups are there, along with a variety of tables with various demographic data. They are not prejoined, but whatever, I can still join them by ID. But the field names in the tables are some kind of census code, e.g. B01001e2.

There is a table that gives a description of what these mean, but how do I translate that into field names for the other tables?

There are relationship classes for the tables, but they don't appear to be actually set.


texts:
I have geotiff file, please see gdalinfo

Size is 28249, 24333
Coordinate System is:
PROJCS["WGS 84 / UTM zone 38N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",45],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32638"]]
Origin = (353455.689043999998830,4691159.864544000476599)
Pixel Size = (0.600000000000001,-0.599999999999992)
Metadata:
  AREA_OR_POINT=Area
  DataType=Generic
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  353455.689, 4691159.865) ( 43d13'13.84"E, 42d21'32.06"N)
Lower Left  (  353455.689, 4676560.065) ( 43d13'27.15"E, 42d13'38.93"N)
Upper Right (  370405.089, 4691159.865) ( 43d25'34.49"E, 42d21'42.89"N)
Lower Right (  370405.089, 4676560.065) ( 43d25'46.26"E, 42d13'49.71"N)
Center      (  361930.389, 4683859.965) ( 43d19'30.43"E, 42d17'41.06"N)
Band 1 Block=28249x1 Type=Float32, ColorInterp=Gray
  NoData Value=256
  Metadata:
    RepresentationType=ATHEMATIC
Band 2 Block=28249x1 Type=Float32, ColorInterp=Undefined
  NoData Value=256
  Metadata:
    RepresentationType=ATHEMATIC
Band 3 Block=28249x1 Type=Float32, ColorInterp=Undefined
  NoData Value=256
  Metadata:
    RepresentationType=ATHEMATIC

My two step to create tiles
1. gdal_translate -ot byte -of vrt source.tiff destination.vrt
2. gdal2tiles -p raster -w leaflet destination.vrt outfupfolder

 image in tiff

 image after gdal2tiles

How to avoid losing color?
please help my

UPDATE_1
So, pct2rgb write error

ct_size = ct.GetCount()
 AttributeError: 'NoneType' object has no attribute 'GetCount'

I did not find an explanation for this error I'm stumped, who can help me


texts: Open Source equivalent to Collector for ArcGIS to use with QGIS/PostGIS?

texts:
From my reading of the Esri web pages the Collector for ArcGIS is tightly tied to that platform.  I would be delighted to learn that I am wrong about this.

I have a strong interest in generic mobile data collection apps for iOS and Android devices with GPSs.  Particularly ones that can collect data while out of coverage, and then sync when they get back in to coverage.  Ideally, they would also be able to cache portions of the data locally for use when off line.

What I envision is a web based application that allows one to specify your data layout, constraints, etc. and which will then generate the the backend database and allow access via a restful API.

What I need to do is enter data from a form on an mobile device, preferably both iOS and android, which will update tables in a backend database, potentially something with GIS extensions such as PostGIS.  That would enable me to easily incorporate the data into QGIS.  It must be able to work offline and many of the areas we work have not 3G coverage.  Ideally you would be able to also cache selected data on the device for offline use.  Map integrations is not important.

I have found some good cloud based solutions but I really want something where I have control of the back end data collection as some of the data is potentially sensitive.  e.g. reporting sighting of rare reptiles which are a smuggling threat.

What alternatives are the for mobile apps to get data into QGIS.


texts: PostgreSQL add outcome of query as column to table "ERROR: more than one row returned by a subquery used as an expression"

texts:
I have a table with boat positions. I am trying to find the the closest waterway to these data points and add the name of this waterway in a new column. I have tried to realize this with the following query:
ALTER TABLE boat_positions
ADD COLUMN wayname VARCHAR;

UPDATE boat_positions SET wayname = (SELECT DISTINCT ON (a) b.name
FROM public."boat_positions" a, public."waterways" b
WHERE b.geom && ST_Expand(a.geom, 100)
ORDER BY a, ST_Distance(b.geom, a.geom) ASC);

With this I run into the error:
ERROR:  more than one row returned by a subquery used as an expression
SQL state: 21000

I understand why this is but I don't have the SQL skills to solve it in a good way. I came up with an indirect solution by first creating a new table and then combining it with the original, as below:
CREATE TABLE wayname AS
SELECT DISTINCT ON (a) a.index, b.name
FROM public."boat_positions" a, public."waterways" b
WHERE b.geom && ST_Expand(a.geom, 100)
ORDER BY a, ST_Distance(b.geom, a.geom) ASC;

ALTER TABLE boat_positions
ADD COLUMN wayname VARCHAR;

UPDATE boat_positions AS a
SET wayname = (SELECT w.name FROM wayname AS b
WHERE a.index=b.index)

However, this solution takes time and memory. Is there better way to do this directly?


texts:
I have a shapefile of points (samples collection locations) with their longitude and latitude and the number of samples collected in each location.
My goal is to create a raster of the number of samples collected per zip-code.
Unfortunately my shapefile doesn't include zip code information.

Is there a way to get the zip code for every sample collection location? I was also thinking maybe there is a way to have that shapefile on top of a US zipped basemap and combine both?


texts:
I am currently trying to serve a Raster Datasheet made from several Geo(TIFF) files but I am unsure about how to go about it as I'm quite new to GIS! I need to host it online so that clients can request tiles from the server in the format: http:///tile/{z}/{x}/{y}/

For example:
http://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/tile/0/0/0 - something similar to this; even better with the REST API.

I am unable to host it on ArcGIS Online as I am not part of an organisation so it will be expensive for a personal uses.

Are there any free methods of hosting raster datasheets as tile servers?

If I create a tile cache using GeoServer and GeoWebCache how can I publish it so online users can navigate to the map and request tiles?


texts: Calculate MBR envelope for different geometry shapes in different rows

texts:
I have some different kinds of geometry (Polygons and line strings) stored in my mysql database. On different rows.

I know I can retrieve the envelope of one row with: SELECT Envelope(zoneShape) FROM geometryTable WHERE refid=4

If there are three different shapes associated with reference id 4, i'm getting three results, obviously.

How can I (preferably with one sql query) get the envelope/MBR of all the geometry with e.g. refid=4?


texts: Calculate perimeter of a large polygon using geodesic calculations

texts: How do I calculate perimeter of a large vector polygon (polygon set) which spans almost half the globe? I guess regular calculations do not produce accurate results in this case.


texts:
I'm trying to add Hawths Tools to by Toolbar selection.  I'm running ArcMap/ArcInfo v.9.3.1. Hawths Tools is not supported in v.10.  I've dowloaded the zip file, extracted the files, adn ran teh .exe file.  When I go to "customize" Tools, add from file, and select under Program files "HawthsTools3.dll", an "Added Objects..." menu pops up that reads: No new objects added.

If you then open ArcMap, select View, then Toolbars, Hawths Tools does not show up.

Can anyone provide assistance?


texts: PyQT5 - How to use QCompleter, QGIS Plugin

texts:
I'm trying to use QCompleter to get something like in the image. The code seems pretty simple, but in my dialog nothing happens.

Is there some problem in this code?

The dialog is created using Qt Designer.

import sys
from PyQt5.QtWidgets import QAction, QFileDialog, QApplication, QWidget, QVBoxLayout, QLineEdit, QCompleter

#

def run(self):
  if self.first_start == True:
    self.first_start = False
    self.dlg = SelectFeatureDialog()

    strList = ['Python', 'PyQt5', 'Qt', 'QML']
    completer = QCompleter(strList)
    self.dlg.lineEdit.setCompleter(completer)

  self.dlg.show()

Are there some errors in these settings?

I'm very stupid, the completer works but I was trying with the lowercase, while in my string's list the first letter is always uppercase.

Now I'm wondering how to ignore the case?