texts: Using qgis python to create a PDF output with print composer to only show the attribute table

texts:
I'm attempting to show the information of a given vector layer that I have by creating a PDF file that shows the attribute table of that same vector layer.

The problem I have is that currently I get an output, but only up to the sixth feature of the attribute table. However I want to display the entire attribute table no matter how many features there are.

I have tried using different layers, same result occurs
I have looked online for a solution extensively. The code that I am using is based of the sources that I have found by looking online.
I am not looking to include any other information other than the attribute table as of now
It is not the number of features or attributes. I've tested a layer with 10 features which should easily fit on the page but no features past the sixth feature are shown.

I believe that something with my code needs to be altered, however no matter how I tinker with it, I can't seem to get more than 6 rows to appear on the PDF file. I'm not entirely sure what needs to be changed as it's always 6, and only 6 rows.

Code

from qgis.gui import *

huff_model = iface.activeLayer()

mapRenderer = iface.mapCanvas().mapRenderer()
c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)
x, y = 0, 0
w, h = c.paperWidth(), c.paperHeight()
composerMap = QgsComposerMap(c, x, y, w, h)

table = QgsComposerAttributeTable(c)
table.setComposerMap(composerMap)
table.setVectorLayer(huff_model)
c.addItem(table)

printer = QPrinter()
printer.setOutputFormat(QPrinter.PdfFormat)
printer.setOutputFileName("C:/Users/Oski/Py_internship/Huff_Model_Oski/Outputs/out.pdf")
printer.setPaperSize(QSizeF(c.paperWidth(), c.paperHeight()), QPrinter.Millimeter)
printer.setFullPage(True)

pdfPainter = QPainter(printer)
paperRectMM = printer.pageRect(QPrinter.Millimeter)
paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
c.render(pdfPainter, paperRectPixel, paperRectMM)
pdfPainter.end()

This code can be run straight from the python console (I used it like this for testing), provided that there is a valid vector layer with an attribute table.

PDF Output

My apologizes for the poor image, it was difficult to capture the entire page with decent resolution. Nevertheless, it's evident what my output looks like; there are only 6 rows shown despite the fact that the vector layer has over 1000 features. I am aware that not all the attributes fit on the page, I'm only concerned with why the remaining features do not appear.


texts: Getting Polyline (Linestring) coordinates from line layer in QGIS

texts:
I have a road network (line layer). I am trying to get list of coordinates (Lat,Long) from road network in the form of geometry LINESTRING as shown below. However, I am only getting two coordinates as shown in second figure.

What I have tried:
lyr = iface.activeLayer()
fl = ['lat','lon'] #Add these fields before executing code

with edit(lyr):
    for f in lyr.getFeatures():
        p = f.geometry().centroid().asPoint()
        _=f.setAttribute(f.fields().indexFromName(fl[0]), p.y())
        _=f.setAttribute(f.fields().indexFromName(fl[1]), p.x())
        _=lyr.updateFeature(f)


texts: Left join resulting in weird output in QGIS's Virtual Layer

texts:
I want to join two layers, a polygon layer "plan" and a point layer "education point". My goal is to join the data from "education point" layer to the "plan" layer if the point is inside a polygon.
Since some polygons don't intersect with any point, and I still want to keep those polygons visible, I used left join. But the result only shows polygons intersecting with a point.
Here's my query:
SELECT pl.geometry,ed.*
FROM "plan" AS pl
LEFT JOIN "education point" AS 'ed' ON st_intersects(pl.geometry,ed.geometry)
WHERE pl."function" in ('primary school', 'middle school')
AND ed."studyArea"='yes'


texts: Training samples for deep learning in ArcGIS Pro 2.5

texts:
I am performing end-to-end deep learning in ArcGIS Pro 2.5, following a similar procedure to the one set out in this tutorial: https://learn.arcgis.com/en/projects/use-deep-learning-to-assess-palm-tree-health/
I am trying to train the model to recognise bats in thermal orthomosaics taken from a drone. The bats are much warmer than their background, so are detectable in the images as bright circular objects with a diameter of approximately 5 pixels, against a dark background (see image 1).
So far, I have been defining training samples by drawing a circle from the center of each bat to just outside the edge of the animal.
However, this is not giving me great results.
Is anyone able to offer advice on the best way to define a bat in my training images?
Should I exclude all of the dark area surrounding the bat?
Should i use the square training sample drawing tool instead?


texts:
I have the following data base with information about rainfall for a set of meteorological stations:

  year  station latitude longitude depto municipio       date rainfall
  1990 11010010   0527 N    7632 W CHOCO     LLORO 1990-01-01     13.2
  1990 11010010   0527 N    7632 W CHOCO     LLORO 1990-01-02     83.1
  1990 11010010   0527 N    7632 W CHOCO     LLORO 1990-01-03     10.5
  1990 11010010   0527 N    7632 W CHOCO     LLORO 1990-01-04      9.1
  1990 11010010   0527 N    7632 W CHOCO     LLORO 1990-01-05      0.0
  1990 11010010   0527 N    7632 W CHOCO     LLORO 1990-01-06     10.4

I am trying to import this to QGIS using the "New layer from delimited text" tool, but QGIS is unable to detect the longitude and latitude even specifying the X and Y coordinates. I know these longitudes and latitudes are not the common WGS84 ones, there is a way to reproject both coordinates?

(I also use R but rgdal commands do not work with the Spatialframe).


texts:
I have two shapefiles, one is polyline, the other polygon. For each polygon I would like to clip the polyline vectors inside. Using arcpy I already found part of a possible solution: https://gis.stackexchange.com/a/19934/6458

import arcpy

arcpy.env.workspace = r'D:\Projects\GDBs\slowbutter.gdb\IPAS'
rows = arcpy.SearchCursor('HspAOI')
for row in rows:
    feat = row.Shape
    arcpy.Clip_analysis('HspWBD_HU12', feat, 'HspWBD_HU12_' + str(row.getValue('NAME')), '')

This however would create a new shapefile per polygon. I have 8000 polygons and would have to merge these using arcpy.Merge_management().

Is there maybe a way to not have to merge 8000 separate shapefiles? I work with arcpy and ogr so the latter might also be an option?


texts: Custom ArcObjects Toolbar Unavailable for Some Users

texts:
I'm responsible for maintaining an extension to ArcMap that consists of a toolbar with a bunch of menus & commands & tools.  It uses the old-style COM registration for all of its UI. Recent installations on our test machines have intermittently refused to put the toolbar in the list of toolbars that users can select.  It is listed for me (a local admin) always; for others sometimes, and for still others never.

Even when the toolbar is not listed, most of the custom commands and such (likely all of them, but there are several dozen and I'm too lazy to take a definitive tally) are available in ArcMap's Customize window.

At first I assumed that the COM registration was failing specifically for the toolbar, but there is a registry entry for HKEY_CLASSES_ROOT\CLSID\MyToolbars.ProgID.MatchesMyToolbars.NetNamespace so I think everything looks okay there.

Where else might I look to find a cause of this issue?

Updates & Fresh Breadcrumbs

Based on Kirk's advice in the comments I confirmed that the toolbar is registered using the categories.exe utility.
A coworker discovered that if he enables the extension that comes in the same library as the toolbar, closes ArcMap, then reopens ArcMap, the toolbar appears in the list.  Deleting the user's profile indeed resets the problem.
A well-meaning sysadmin found that running CCleaner clears up the problem, but they ran it wholesale and made no log of what was changed.  On the one remaining test machine that wasn't CCleaned, I re-did that cleaning incrementally to determine which piece of flotsam causes the problem but came up empty-handed.
All but one of the test machines (the one that I cleaned manually in (3)) have been rebuilt from a new image and the problem persists.  If anything its intermittence has decreased, and the toolbar is more likely not to be available.
After the rebuild I was able to determine with a moderate level of confidence (still some intermittence) that the installation protocol I use on our test computers contains bugs.  Working on a potential fix now.


texts: Esri JavaScript ArcGIS API - Developing a FeatureMap Attributes Filter Widget

texts:
I'm currently trying to develop a FeatureMap Attributes Filter for Esri basemap with a featuremap on top of it, but I am struggling with some of the new stuff as I am fairly new to this type of programming (Python, data arrays analyst here).

My problem and question:

How do I input my current FeatureMap so the widget knows what to analyse?

Here is my widget code:

/////////////////// Widget
@subclass("FilterWidget")
class FilterWidget extends declared(Widget) {
  constructor(params?: any) {
    super();
  }

  // Variable input for the filtering
  @property()
  @renderable()
  FirstAttribute: string;

  // Variable input for the filtering
  @property()
  @renderable()
  SecondAttribute: string;

  // Variable input for the filtering
  @property()
  @renderable()
  ThirdAttribute: string;

   define the property dependencies
  @property()
   readOnly: true;
   dependsOn: ["FirstAttribute", "SecondAttribute", "ThirdAttribute"];

  //// Drop-down menu layout
  render() {
    var sqlExpressions = ["ID > 0", this.FirstAttribute, this.SecondAttribute, this.ThirdAttribute];

       var selectFilter = document.createElement("select");
       selectFilter.setAttribute("class", "esri-widget esri-select");
       selectFilter.setAttribute("style", "width: 275px; font-family: Avenir Next W00; font-size: 1em;");

       sqlExpressions.forEach(function(sql){
         var option = document.createElement("option");
         option.value = sql;
         option.innerHTML = sql;
         selectFilter.appendChild(option);
       });
       //view.ui.add(selectFilter, "top-right");

  //// Filtering computation
       function setFeatureLayerViewFilter(expression: any) {
           view.whenLayerView(layer).then(function(featureLayerView: any) {
             featureLayerView.filter = {
               where: expression
             };
           });
         }

       selectFilter.addEventListener('change', function(event) {
       setFeatureLayerViewFilter((event.target as HTMLInputElement).value);
     });

     return (
      <div class="browser esri-widget esri-widget--panel">
      </div>
    );
  }
}
export = FilterWidget;

and here is the following I write in my main.ts file to call the widget

var filterByEjer = new FilterWidget({
  FirstAttribute: "Ejer = 'Ejer'",
  SecondAttribute: "Ejer = 'Uoplyst'",
  ThirdAttribute: "Ejer = 'Kommunal'"
});

view.ui.add(filterByEjer, "top-right");


texts: Returning Geometry - ArcGIS REST service

texts:
Query Esri service won't return geometry.
I'm using a query on ArcGIS REST service, the first link returns geometry, but with the same query parameters the second link won't return geometry.

First link with geometry
the second link without geometry, but it should be with geometry.

They both are using the same parameters, and I can assure you that the second link features have geometry (for instance).
So,
how to return geometry on the second link?


texts: Getting slope from DEM.tif using Spatial Analyst

texts: I have a DEM.tif file and I want to get the slope. I know I use the Slope tool from Spatial Analyst. But, I'm confused which output_measurement to choose from (i.e degrees or percent rise). Basically, I am trying to find areas where the slope < 5%.


texts: 'VCRUNTIME140.dll and MSVCP140.dll' missing when loading QGIS

texts:
I have downloaded QGIS 2.18 on Windows 10 without any problem but and when I start the application, I am getting the below error:

  missing VCRUNTIME140.dll and MSVCP140.dll


texts: netcdf daily data to monthly raster layers

texts:
I have a netcdf (.nc) file of 95GB which consisted of daily data for a 50 years period. I tried to process this with ArcGIS (Multidimension tools – MakeNetCDFRasterLayer command) and it was not functioned. Then I tried with QGIS by adding .nc file as a raster layer and ‘save as’ as Gtiff. This also crashed at the middle of the process.
.nc file data are in daily basis and I want to create separate monthly raster layers by using daily data.
Can someone help me solve this?


texts:
I just computed multiple ring buffers for a point in EPSG 4326 with a 100m interval for a maximum distance of 1000m. I was observing a calculation time of approximately:

50s within ArcMap;
4min in ArcGIS Pro; and
40s as a processing service on ArcGIS Server 10.5.1.

I would like to see the same performance in ArcGIS Pro as well as ArcMAp because I am shipping a toolbox with this arcpy command:

def createDistances(incr,maximum):
    #as we have the raster now, we will need a multiple ring buffer:
    distances = []
    for dist in range(incr,maximum+1,incr):
        distances.append(str(dist))
        distancesString = ";".join(distances)
    return distancesString

distancesString = createDistances(100,1000)
arcpy.analysis.MultipleRingBuffer("Your single point feature", r"your ring buffer name", distancesString, "Meters", "distance", "ALL", "FULL")

Are there ways to increase calculation times in ArcGIS Pro?

I am running this on a X270 with this spec:

Intel(R) Core(TM) i7-7600U CPU @ 2.80GHz, 2904 MHz, 2 Cores, 4 logical processors
16Gb Ram, SSD

no dedicated GPU.


texts:
I have compared some distance calculations between geopy and PostGIS ST_Distance_Sphere, and it seems that PostGIS ST_Distance uses a calculation that is similar to the great circle distance rather than Vincenty.

Geopy Vincenty (Most accurate) using geopy

>>> from geopy.distance import great_circle, vincenty
>>> p1=(32.083659,34.806480)
>>> p2=(32.0843373,34.8065365)
>>> vincenty_distance = vincenty(p1, p2).meters
>>> print vincenty_distance
75.4044140365

Great Circle (less accurate) using geopy

>>> great_circle_distance = great_circle(p1, p2).meters
>>> print great_circle_distance
75.6324224394

PostGIS ST_Distance_Sphere, Unknown Algorithm:

=> SELECT ST_Distance_Sphere(
          ST_SetSRID(ST_MakePoint(34.806480, 32.083659), 4326),
          ST_SetSRID(ST_MakePoint(34.8065365, 32.0843373), 4326));
 st_distance_sphere
--------------------
       75.611223453
(1 row)

Does PostGIS have a Vincenty Distance calculation?


texts: Elevation graph shows all the profiles starting at the same point eventhough one of them starts farther inland?

texts:
I am trying to compare 3 interpolated lines to analyze sediment movement after beach nourishment. Two of the lines start at the same point, but a the 3rd one starts farther inland. However, when I put them all in a graph they are all shown starting at the same point.

Is there a way I can fix this problem?

I tried editing the graph properties and it doesn/t work. All shapefiles and DEM are in the same TX state plane projection. I can't think of any other reason why they would plot correctly. [X axis distance/Y axis elevation, both in ft. Pink and yellow lines in the graphs should start at about 125 ft.


texts: Implicitly convert Geometry to WKT string in PostGIS

texts:
Suppose I have a table defined as follows:

CREATE TABLE gtest (name varchar, geom geometry);

To insert, I am able to simply do this:

INSERT INTO gtest
VALUES (
    'Polygon',
    'SRID=4326;POLYGON((0 0,1 0,1 1,0 1,0 0))'
);

I don't have to wrap the WKT string in the function ST_GeomFromText() because PostGIS has an implicit cast that does so. This is nicely explained by @JGH here

By using the PostgreSQL command \dC, the defined casts can be listed, including:

                                     List of casts
     Source type         |         Target type         |      Function      |   Implicit?
-------------------------+-----------------------------+--------------------+---------------
text                     | geometry                    | geometry           | yes
geometry                 | text                        | text               | yes

I'd like to make it, so I can simply do a SELECT * FROM gtest, and have the results of the geometry column be implicitly converted into WKT. Currently, it will just display them as WKB.

First, I tried creating a new cast as follows:

CREATE CAST (geometry AS text) WITH FUNCTION st_astext(geometry) AS IMPLICIT;

This returned an error, since a cast from geometry to text already exists (as seen in the table).

I then tried ALTER EXTENSION postgis DROP CAST (geometry as text); and then DROP CAST (geometry as text);, and was then able to create the new cast:

                                     List of casts
     Source type         |         Target type         |      Function      |   Implicit?
-------------------------+-----------------------------+--------------------+---------------
geometry                 | text                        | st_astext          | yes

This still didn't work however, as when I do a select, I still get the results in WKB.

Firstly, is this possible? Am I just doing something wrong?
Secondly, would any geometry functions break by adding this implicit cast?


texts: Setting up lat/long graticule with 1 min output in QGIS?

texts:
I would like to set up a lat/long graticule in QGIS with output being a line shapefile with values in Degrees and Minutes, and grid lines shown at every even 1 minute (e.g. 13°23' W, 13°24' W, 13°25' W  ..., and likewise for the other axis), and if possible have these values written in this exact format in the attribute table of the output file.

It seems that there are more ways to how to set up the graticule, e.g. in the program itself (Vector -> Research Tools -> Vector Grid, or via MMQGIS plugin) or in the Print Composer, however none of these methods provide a straight forward output as desired; only in decimal numbers.

How can I do it?

QGIS version 2.18.14 and my project and its layers are in EPSG:3857.


texts:
I have written a script which takes two inputs, intersects them, creates a new field in the output feature class and joins it back to one of the input feature classes. I would then like to select the records in the joined 'Intersect' field that contain the string 'Yes'. Do I need to convert the input from arcpy.GetParameterAsText(2) to a string and use this in the expression? My code is below.

import arcpy

arcpy.env.overwriteOutput = True

inFeatures = (arcpy.GetParameterAsText(0), arcpy.GetParameterAsText(1))
outFeatures = arcpy.GetParameterAsText(2)
joinField = "Ref"

arcpy.Intersect_analysis(inFeatures, outFeatures)

arcpy.AddField_management(outFeatures, "Intersect", "TEXT", field_length =
50)

arcpy.CalculateField_management(outFeatures, "Intersect" , '"""Yes"""',
"PYTHON")

arcpy.AddJoin_management(arcpy.GetParameterAsText(0), joinField, outFeatures,
joinField)

arcpy.SelectLayerByAttribute_management(arcpy.GetParameterAsText(0),
"NEW_SELECTION", "\".[Intersect]\" = 'Yes'")


texts:
What is the problem with gdal info for sentinel 2 level 2A?

gdalinfo SENTINEL2_L2A:S2A_MSIL2A_20181227T125301_N0211_R052_T24LXN_20181227T145100.SAFE/MTD_MSIL2A.xml
ERROR 1: Invalid syntax for SENTINEL2_L2A:
gdalinfo failed - unable to open 'SENTINEL2_L2A:S2A_MSIL2A_20181227T125301_N0211_R052_T24LXN_20181227T145100.SAFE/MTD_MSIL2A.xml'.


texts: Raster to Polygon tool results in reduced area

texts:
I have a large raster dataset within a file geodatabase (FGDBR format raster) which I have converted to polygon using the 'Raster to Polygon' tool.

Cell size of original raster = 25m.

Cell size of polygon = 25m.

Total area of original raster = 10,423,764 ha.

Total Area of polygon dataset = 10,235,534 ha.

Total area of raster when values extracted then added together (see 5th paragraph)                    = 10,235,534 ha.

Projected coordinate system = British National Grid

The area of the resulting polygon dataset is 1.8% less than area the cells of the original raster. However, visually; there is no difference between the two datasets (checked by overlaying them in ArcMap), except when 'zoomed' right out, then you can see slithers of the original raster cells outside the polygons. When scale is less than 1:200,000, the slithers disappear.

I calculated the original raster area by multiplying the cell count by the area of each square. I did not check the 'Simplify Polygons' option when using the tool.

I've also created the polygon dataset in different ways, such as using the 'Extract by Attributes' tool to extract separate parts of the raster and converting each part to polygon. But this also results in an area of 1.8% less, but no visual difference except at 1:200,000 scale or more.

I've also used Raster Calculator to subtract the extracted parts of the raster from the original raster, but this results on a raster with a value of 0 throughout.

Can anyone explain what is going on here please?


texts: What are the top best webmapping API clients?

texts:
Possible Duplicate:
  Comparison of JavaScript mapping libraries?

I am looking for the top APIs that I can use to display maps from a number of map servers/sources. (like openlayers and openfscales)
I already know OpenSlales and Openlayers. So are there any other popular APIs for other languages (C++,Java,... ) ?


texts:
Is it possible to distinguish (with ArcPy) a featureclass in an RDBMS (SQL Server, in this case) that is a table vs a registered spatial view that was created with the command line tools?

Spatial view is stored in an ArcSDE 10.0 geodatabase, and arcpy version <= 10.3.1.


texts:
I can easily open and type all the appropriate things in the various tools, but when I execute the tools, it just shows an hourglass in the results window and nothing never happens. This is the same for all tools.
Any idea how to fix it?


texts: Configuring spatialite database access for Python 2.6.6

texts:
I'm working on a GIS project, and I would like to implement and test some geo-spatial algorithms in Python. For this purpose, I will not only need sqlite, but also spatialite, in order to store and query the location data. Now I've tried to install the pyspatialite package, but no matter what Python version I tried (I tried all versions from 2.6 to 3.3), the pip keeps insisting, that none of the existing pyspatialite packages are compatible with my version of Python.

If I try to do this using easy_install, I get a traceback and an error:

AttributeError: MSVCCompiler instance has no attribute 'compiler'

And that also occurs, if I try to install the package manually, by executing the setup.py file.

From what I've already searched, some people suggest to connect to a spatialite database somehow using sqlite and loading extension, but frankly I have no idea how to do it, and couldn't understand any of these answers. I will be really grateful if someone here is able to propose a solution in a clear, step-by step way, as I'm not a very experienced Python programmer yet. Thanks in advance.

UPDATE:

Another attempt, this time with Python 3.3.5. The following code:

import sqlite3
conn = sqlite3.connect(":memory:")
conn.enable_load_extension(True)
conn.execute('SELECT load_extension("libspatialite-2.dll")')

yields:

Traceback (most recent call last):
  File "<pyshell#10>", line 1, in <module>
    conn.execute('SELECT load_extension("libspatialite-2.dll")')
sqlite3.OperationalError: %1 is not a valid Win32 application.

And again, I don't seem to resolve this on my own. I'm really stuck here, and slowly running out of ideas. If you can make it work on any Python version, I'd be REALLY relieved.

#

UPDATE #2

Sorry, the error above was due to something else, I reinstalled Python and pysqlite, and we are back with the old error. There are two options now:

1) I compile using from pysqlite2 import dbapi2 as sqlite3. In this case code is the following:

from pysqlite2 import dbapi2 as sqlite3
conn = sqlite3.connect(":memory:")
conn.enable_load_extension(True)
conn.execute('SELECT load_extension("DLLs\\libspatialite-4.dll")')
curs = conn.cursor()

In that case, the error is:

Traceback (most recent call last):
  File "C:\Users\mszydlowski\Desktop\Project\sqlite.py", line 3, in <module>
    conn.enable_load_extension(True)
AttributeError: 'pysqlite2.dbapi2.Connection' object has no attribute  'enable_load_extension'

2) I compile using import sqlite3. In this case code is the following:

import sqlite3
conn = sqlite3.connect(":memory:")
conn.enable_load_extension(True)
conn.execute('SELECT load_extension("DLLs\\libspatialite-4.dll")')
curs = conn.cursor()

And the error:

Traceback (most recent call last):
  File "C:\Users\mszydlowski\Desktop\Project\sqlite.py", line 4, in <module>
    conn.execute('SELECT load_extension("DLLs\\libspatialite-4.dll")')
OperationalError: The specified module could not be found.

even though the file is certainly in there, and I did use the double backslashes, like I was already advised.

Any ideas?


texts: Download layers from arcGIS online in python

texts:
How can I get layers/maps from arcgis online in Python? I mean to download layers from Arcgis online in python, just like publishing.
I've searched about it in the documentation of Arcgis 10.2.2, but I  only found publishing.

At the end, I want to write a Python script that will download, automaticly, layers\feature classes from the arcgis online, by specific areas.
Please help me.


texts:
Possible Duplicate:
  How to convert .shp into .gpx using QGIS?

I'm trying to convert a shapefile with some random points into a gpx file, but I'm having some problems: when I upload the gpx file into my Garmin Etrex, the points doesn't have coordinates. Why?

ps: I'm working with the EPSG 3003, and I've regulated the Garmin too.

ps: I'm using QGIS 1.8 Lisboa on Mac OsX Mountain Lion. The main problem is that the GPS Tools is not working properly


texts: Large number of distribution maps - Batch maps / Atlas

texts:
I try to make distribution maps for about 500 species.
In the orginal file all species / points are in one shapefile, via the tombio productivity tool-pluging I separated the species and created simple maps.
(https://www.youtube.com/watch?v=Iq6UWPcH46k&feature=youtu.be)

However I need the maps with title, number of features etc. in it and publication quality.

I found
Batch maps (atlas plugin style) in same location
and
How to create an Atlas based on point attributes?

but none of these seem suitable for the high number of species.

The result should be the ca. 500 maps with the possibilty to add text fields with expressions from the shapefile.

Is there "simple" solution in Qgis via the atlas generation?


texts: Multiplying values within ArcPy update cursor?

texts:
This code is meant to apply an update cursor to a field called "New_Sq_Mile" when I am populating the field with the area of the polygons. Unfortunately the units are in meters and I need to do a calculation to change those values to miles before they are updated into the table. Here is what I have so far.

cursor = arcpy.UpdateCursor(FinalErase, ["SHAPE@AREA", "New_Sq_Mile"])

for row in cursor:
    row[1] = row[0]
    update.row

del row
del cursor

Is there a way I can do a calculation in the for loop such as

row[1] = row[0]/1606.344


texts:
I have raster tiles in mbtiles format.
I also have vector tiles in GeoJson format.
My goal is to combine both tiles so they will overlap each other, the vector tiles will be on top of the raster tiles.
I am using QGIS to import both layers of tiles and can't find a way to make them overlap.


texts: Using field mappings to not carry over fields

texts:
I've tried multiple times to understand field mappings and continue to be rather baffled. I assume it can be used to perform a geoprocess (in my case it will be a spatial join) and limit the fields in the output feature class. I'd like have the output only contain the fields generated by the spatial join (Join_Count, etc). I've tried creating blank field mappings objects and using it, but that's no good. I'd rather not take the extra step to delete all fields, as in large data set it will take some time.

Can anyone provide help?


texts: Styling based on hex color codes within table?

texts:
I have a table of polygons named 'testtable'.  Within it is a column 'color', which has hex codes.
I want to style the polygons based on the hexcode values.
Is this possible?
I was trying something like this, but it's obviously not working..
# testtable{
  polygon-fill: {color};
}


texts: Is there a measurement standard e.g. for coast lines?

texts:
When watching discovery & documentary shows, coastlines are apparently defined to be just long enough to fit some explanation or a discussed topic.

The coastline paradox is a well-established fact describing how the measurement is defined by the precision of measurement.

By specifying some sensible alignment, precision and assumptions for e.g. tides, one could imagine a well-defined algorithm for proper & comparable measurement (e.g. by rasterization).

So is there any established standard for measuring coast lines?


texts: PyQGIS using selected layer from a combobox

texts:
With PyQGIS 3, I defined a comboBox to list layers from a QGIS project. I want to use the selected one to execute a treatment. I tried comboBox.currentData() class but it doesn't seem to be the good way.

How can I handle the selected layer from the comboBox as a vector layer ?

Maybe with comboBox.itemData() or comboBox.currentIndex() classes ?

# ComboBox
layers = QgsProject.instance().layerTreeRoot().children()
self.dlg.comboBox.addItems([layer.name() for layer in layers])
# Not the good way
selectedLayer = self.dlg.comboBox.currentData()

# Vector Layer
mySelectedLayer = QgsVectorLayer(selectedLayer, 'polygon', 'ogr')


texts: Merge all features in polygon shapefile with gdal

texts:
How do I apply using this same concept,  but simply merging all features into one?

ogr2ogr output.shp input.shp -dialect sqlite -sql
"SELECT ST_Union(geometry), dissolve_field FROM input GROUP BY dissolve_field"


texts: Finding widest portion of river bank using ArcGIS Desktop?

texts:
I am having a long portion of river. Here I am displaying part of it. In this polygon which is transparent black.

How do I find the wider track (bank to bank perpendicular from the center line) using ArcGIS Desktop?


texts: I am using QGIS 2.18 to identify points that fall out of a specific boundary but the identify features is not working. When I right click, it brings me the polygon feature. So i turned off the polygon layer and highlighted the layer i am working with as described in QGIS identify features button issue but it is still not working. I found the identify mode in the map tools but cant seem to figure out what: "You probably want the Top down, stop at first" means. It is either it does not bring up the table at all or the identify results table is blank. I have also tried this in QGIS 3.4 and got the same results. Some points are identified while some are not.


texts:
File geodatabase (fgdb) includes numerous file geodatabase tables. As far as I know they exist as dbf files, but are within a Database.gdb.

In ArcCatalog, the file pathway resembles C:\Users\...\Database.gdb\Stats_AA.

How to read all of these dbf files into R (a statistical software)? What is the correct pathway to supply? The function used is read.dbf (in the foreign package).

Variants of

test<-read.dbf(file="C:/Users..Database.gdb/Stats_AA.dbf")

and

test<-read.dbf(file="C:/Users..Database/Stats_AA.dbf")

don't work. What's the correct "form" of the file name to be used, or, do I need to export all of the file geodatabase tables into some other form or location?


texts: Adding WMS into the Google Maps or Mapbox

texts: Is it possible for me to get the data from https://firms.modaps.eosdis.nasa.gov/web-services/ and add it as a layer in the maps so that whatever data I am getting from the server can be directly shown on my web page using Mapbox or Google API. I am new to this and like to integrate it in my project.


texts: Converting JSON to GeoJSON MultiPoint using `ogr.CreateGeometryFromJSON`

texts:
I am trying to convert a JSON array to a geom in Python, by first converting it to a GeoJSON.

I have point coordinates in a JSON array,

[{'lat':100.0, 'lng':0.0},{'lat':101.0, 'lng':1.0'}]

that I want to convert to GeoJSON to then use as a geom in Python.

I convert the JSON array using a Python function,

def json2geojson(items):
    import json
    return json.dumps({
            "type": "MultiPoint", "coordinates":[ [ feature['lng'], feature['lat']]
                      for feature in json.loads(items)
            ]})

my_geojson = json2geojson(my_json)

which gives me this:

{ "type": "MultiPoint",
    "coordinates": [ [100.0, 0.0], [101.0, 1.0] ]
    }

Then I try to convert a geometry from the GeoJSON,

from osgeo import ogr

geom = ogr.CreateGeometryFromJson(my_geojson)

which at first glance appears to work

<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at 0x113ea9ea0> >

however something clearly goes wrong

print "%d,%d" % (geom.GetX(), geom.GetY())
0,0
print "Geometry has %i points" % (geom.GetPointCount())
Geometry has 0 points

Does anyone know why this is not working or the best way to do this type of conversion?

Edit:

Code works for a single point:

from osgeo import ogr

geojson = """{"type":"Point","coordinates":[108420.33,753808.59]}"""
point = ogr.CreateGeometryFromJson(geojson)
print "%d,%d" % (point.GetX(), point.GetY())

From http://pcjericks.github.io/py-gdalogr-cookbook/geometry.html#create-geometry-from-geojson


texts: foreground arcpy.mapping.MapDocument("CURRENT") in python script in single ArcMap produces createobject errors

texts:
I'm running ArcGIS 10.3.1. and desire to list layers in the current mxd. I'm using the "CURRENT" method...

arcpy.mapping.MapDocument("CURRENT")

...within a python script tool, accessed via a single ArcMap session with the environment setting set to foreground processing. The usage is correct (identical to the help) and has worked for me last year. why would a foreground process fail to find the single mxd that it's being run from?

Link Esri ArcGIS 10.3 Help: Map Documents

    import arcpy
    mxd = arcpy.mapping.MapDocument("CURRENT")
    for lyr in arcpy.mapping.ListLayers(mxd):
        arcpy.AddMessage(lyr.name)
    del mxd
    arcpy.AddMessage(" ... current mxd ")


texts:
Note the attributes selected.
I would like to select the attributes that cross two polygons.
Note the darker polygons represents a overlayed polygon.
I can do this manually, but is there a tool in ArcGIS Pro that can do this for me?


texts: Transform projection from EPSG:4326 into EPSG:900913

texts:
.

i've got this problem about projection. .i just wanted to know what am i missing. .

i already displayed map base layer using openlayers through

base_layer =  new OpenLayers.Layer.WMS( "OpenLayers WMS",
           "http://vmap0.tiles.osgeo.org/wms/vmap0?", {layers: 'basic'} );

i also displayed my vector layer through

var layer = new OpenLayers.Layer.Vector("Butuan");
    map.addLayer(layer);
    var format = new OpenLayers.Format.GeoJSON();
    var feat = format.read({{my_geojson|safe}});
    layer.addFeatures(feat);

and my position is

map.setCenter(new OpenLayers.LonLat(125.54915811889646, 8.90855351100384),11);

which positions the view as i desired it to be.
when i tried to change my base layer to

base_layer =  new OpenLayers.Layer.OSM("OpenStreetMap");

the center of the map is redirected to another location. .i guess that this problem pertains to the projection. i tried to change the center into

map.setCenter(new OpenLayers.LonLat(125.54915811889646, 8.90855351100384),11).transform(geographic, mercator);

which is defined as

var geographic = new OpenLayers.Projection("EPSG:4326");
       var mercator = new OpenLayers.Projection("EPSG:900913");

no changes has been noticed. .the position of the center is still redirected to another location. .

Can anyone suggest solutions on how to solve this problem?

I am using django framework. . vectorformats and simplejson to convert the retrieved spatial data into geojson from postgis. .


texts: Retrieving Wikipedia's coordinate information?

texts:
I'd like to add icons to my map showing locations which have Wikipedia entries.

The WikiLocation API is a great resource for determining Wiki entries near an XY location, but it has limitations (20km radius, 50 points only) which make it unsuitable for my purposes.

Ideally I need the map to show ALL locations with Wiki entries, so I guess I'll need to download the Wikipedia dump file and parse it for geocoded entries.

Where is this information found in the Wikipedia dump files?

The Wiki database may be downloaded from http://en.wikipedia.org/wiki/Wikipedia:Database_download but I can't see which file(s) I need to search for the coordinates.

For example, the entry for Sydney includes:

Coordinates:    33°51′35.9″S 151°12′40″E

Where is this found in the dump files? Another way to phrase the question is "how can I determine which Wiki entries are geo-coded?"


texts: Cartopy image tiles - determining possible zoom level in advance?

texts:
In Cartopy I'm using an image tile like so:

terrain_tile = cimgt.StamenTerrain()
ax.add_image(tile, zoom)

By default I've set the zoom level to 14 and that works perfectly for my first test case.

However, for some locations I get "A non-empty list of tiles should be provided to merge". If I change the zoom level to 13 it then works fine.

This seems to be location-based. With some locations having tiles available with a higher zoom level.

Is there a way to determine in advance which zoom levels are available for a particular location?


texts:
I have a satellite image. I want to add this image to Autocad, and keep its original coordinates.
How can I import an image to Autocad 2014 and keep its original coordinates?


texts:
I am trying to export selected feature using PyQGIS.
Till now I have tried these codes. After using this I can select feature but don't know how to export selected feature to new shapefile.
canvas = iface.mapCanvas()
AllLayers = canvas.layers()
for i in AllLayers:
    it = i.getFeatures(QgsFeatureRequest().setFilterExpression(u'"Country" = \'India\''))
    i.setSelectedFeatures([f.id() for f in it])
    print "Filter Applied"

Have tried this code: but this is just creating a duplicate of source file (I need selected only)
_writer = QgsVectorFileWriter.writeAsVectorFormat(i, r"C:\Users\XYZ\Desktop\NewFile.shp", "utf-8", None, "ESRI Shapefile")

If anybody knows how can I export shapefile using PyQGIS?


texts: How to convert a String variable to a Raster datatype using Calculate Value

texts: I have a GP service I would like to consume in a Flex application, but one of the input data types can be a raster. ArcGIS Server does not support a raster data type as an input parameter, but according to their documentation here - Input and Output Data Types - you can create a variable in your model/script of type String and convert that into the data type your model requires. In their example they use an Areal Unit, but I need to convert to a raster dataset/raster layer. How would I implement this?


texts: Convert PolygonZ to ordinary Polygon Shapefile

texts:
I'm using QGIS version 2.18.11

When looking in properties > metadata of a layer it says PolygonZ. I would like it converted to Polygon instead so I can merge the PolygonZ shapefile with another. At the moment I generate an error when using vector > data management tools > merge...hence if there is an alternative means to merge the two layers which would be easier than converting that would be fine.


texts: Error generating chart: User memory limit exceeded, when changing dataset from T1_TOA to T1_SR

texts:
So, I got the code to work for T1_TOA band, but then when I switched it over to T1_SR band it gave me this error: User memory limit exceeded.
So, I poked around the forum and found the code to limit my date frames, but I still get the error.
Not quite sure what I'm missing.
// Import the Landsat 8 SR image collecton.

Map.addLayer(study_area, {color: 'Black'}, "study_area");

function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = (1 << 3);
  var cloudsBitMask = (1 << 5);
  // Get the pixel QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
                 .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
  return image.updateMask(mask);
}

var dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterBounds(study_area)
                  .map(maskL8sr)
                  .select('B5' , 'B4');
print(dataset);

var startDate = ee.Date('2013-05-01'); // set analysis start time
var endDate = ee.Date('2013-9-30'); // set analysis end time

// 'Reduce' stack of images by using the median value for each pixel and clip to study area
var output_bands = dataset.reduce(ee.Reducer.median()).clip(study_area);
print(output_bands)
//
var index = output_bands.normalizedDifference(['B5_median', 'B4_median']); // NDVI calculated using NIR and red bands
print(index); // DISPLAY metadata

// DISPLAY output band
Map.addLayer(index, {}, 'NDVI');

print(ui.Chart.image.series({
  imageCollection: index,
  region: study_area,
  reducer: ee.Reducer.median(),
  scale: 30
}).setOptions({title: 'Cloud-masked NDVI over time'}));


texts: Deploying Application With C# GDAL/OGR Bindings

texts:
I developed a relatively simple web app that uses the C# OGR bindings to do a couple of simple geometric functions; buffer, clip...etc.  I developed this application on my local machine, which has FWTools installed, so getting everything to play nicely wasn't very difficult.

Naturally, when I deployed this application to the server nothing worked, undoubtably due to the all of the missing associated dependencies for the C# bindings.  Is there anyway I can include all needed files in the bin folder of this application?  Or is that more trouble than its worth, and I should just install FWTools on the server?


texts:
For my Bachelor's thesis I try to reconstruct a road betweeen two ancient cities.

One method I'd like to try is least cost path analysis with QGIS. I've worked with QGIS before, but mostly only basic stuff.
Can someone explain how least cost path analysis works? I know there are already a lot of threads on this, but I don't really understand.

I downloaded a DEM layer (SRTM satellite data) and transformed it into a slope layer and an elevation layer. Then I used the plugin called "Least Cost Path" on both, but got different results (which isn't suprising I think).

I'd like to take both, slope and elevation, into account for my least cost analysis. How can I combine these two layers? Is there anything else I have to do or any other layers I need?


texts: ST_ShortestLine gives different results for 4326 and 900913 projections

texts:
I'm trying to get my head around how ST_ShortestLine (and ST_Distance and ST_LineLocatePoint) work with different projections. The result for the shortest line differs depending on which projection I use (either 4326 or 900913).

This is the result visualised as 4326 (pink is ST_ShortestLine on geometry left at 4326, blue is ST_ShortestLine on geometry transformed to 900913)

And this is the same result visualised as 900913 (pink is ST_ShortestLine on geometry left at 4326, blue is ST_ShortestLine on geometry transformed to 900913)

As you can see, they are both the shortest line depending on the projection you visualise in. So which one is actually the correct 'shortest line' in reality? The actual thing I'm visualising here is the nearest point on a road to a bus stop – the small difference in output is actually a big deal for geographic accuracy.

If it helps, the final output for this data will be in the 900913 projection (D3.js). I'm happy to go with the transformed 900913 output, I just want to understand what's happening and which one is actually truly shortest.


texts: create a pie shaped object in shapely and export to gis polygon object

texts:
I would like to create a pie chart-shaped polygon with concentric layers like so:

I can create this simple object with:

from math import sin, cos, radians
from shapely.geometry import Point,LineString

#variables for buffers
coordinates = (0,0)
point = Point(coordinates)
buffers=[]
outerRadius=100
numberOfCircles=5
radiusStep= outerRadius/numberOfCircles
innerRadius=radiusStep
pi=3.1415

#create the buffers
while innerRadius<=outerRadius:
    buffers.append(point.buffer(innerRadius,200))
    innerRadius+=radiusStep

#variables for linestrings
centerX = coordinates[0]
centerY = coordinates[1]
coordX=[]
coordY=[]
lines=[]

#create linestrings for circles
for i in xrange(1,361,45):
    xCoord = centerX + outerRadius * cos(radians(i))
    coordX.append(xCoord)
    yCoord = centerY + outerRadius * sin(radians(i))
    coordY.append(yCoord)
    line = LineString([(centerX,centerY),(xCoord,yCoord)])
    lines.append(line)

My question is this:

How can I take multiple shapely objects of different geometries and create a geofile (geojson, shapefile, etc...) so that each section of the "pie" is treated as a separate polygon.


texts:
I made a map with qgis. Unfortunately when I want to render the map with d3 map renderer, it shows weird structures. The map is in "1:116'622'714" scale. Could this be a problem?

The website should display this:

But when I use the Ginzburg VIII, it looks like this:

Same with the Winkel Tripel


texts: QGIS changing the attribute input fields order

texts: I have to draw several stuff and input the data, all the layers have the same data fields but some requires 2\5 attributes to be filled. Is there a way to the change the order of those fields to be on top so I get to input them easier without changing the order of the field themselves?


texts: Copying domain (and only domain) between Esri geodatabases?

texts:
Sometimes I have the need to copy a domain(s) from one geodatabase to another (without the featureclass or table).

In the past, I have accomplished this in a few different ways, but each method is not as straightforward as I think it could be (or maybe there's a CopyDomain tool that I've missed).

What method do you use to accomplish this?


texts: Create Contours From a DEM without 3D or Spatial Analyst

texts:
I'm wondering if there is a way to create contours from a DEM using ArcGIS without 3D or Spatial Analyst Extensions.

I know that there are a few questions out there pertaining to this (How to create DEM from contours without ArcGIS Spatial Analyst licence?), but I'd like to specifically ask if there is a workaround for this. Like I said this is more of a curiosity type question. My organization is looking to cut budget and I personally would like to cut down on some of the ArcGIS licensing cost. Major projects that we do here are elevation related and we need this functionality. Like I said, QGIS can complete it, but some clients request the work be done in ArcGIS.

Currently I'm trying a mix of DEM to points to contours, but with no luck so far. Crunching a DEM to points is a long resource trying task so it's not optimal, but I think it may work. But after that I really need to find a way to go from point to contour. Is there a free ArcGIS add-in with a functionality like this?


texts: Clipping vector by iterating through input clippers using PyQGIS?

texts:
I have a python script that I want to work in QGIS 3.0.

#importing all necessary modules
import processing, os, glob

result_path =("/home/brent/Documents/Hansen_new_data/Regions_Min/Province/TI/CADT/") #specifying the result path/directory
os.chdir("/home/brent/Documents/Hansen_new_data/Mindanao_provinces/") #changing the directory where the input files are located
clipper_path ="/home/brent/Documents/Hansen_new_data/Mindanao_provinces/"

for fname in glob.glob("*.shp"):
    name = list(os.path.splitext(fname)) #convert an iterable to list
    name[0] = name[0] + ""
    newName = "".join(name)
    print ("processing " + newName)
    overlay = newName.replace(".shp", "_clipped.shp") #indicating the ouput shapefile by replacing some of its letter strings
    input_vector = ("/home/brent/Documents/Hansen_new_data/Regions_Min/TI/CADT.shp") #specify the input vector to be clipped
    #processing.runalg("qgis:clip", input_vector, fname, result_path + overlay) #run the clip algorithm
    alg ="qgis:clip"
    params = {"INPUT": input_vector,"OVERLAY": clipper_path+fname,"OUTPUT":result_path + overlay}
    processing.run(alg, params)
print("Finish")

Please take note that the commented line (#processing.runalg("qgis:clip...)
worked in 2.18. I tried to follow the new convention of QGIS 3 processing in this link.

This is the resulting error from the console:

Traceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 22, in <module>
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 105, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 183, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: There were errors executing the algorithm.

Is there something I missed with my python script?


texts: Changing Symbology of Multiple Layers using ArcMap/ArcPy?

texts: Is there a way to change symbology for multiple layers in ArcMap or ArcCatalog. I want to change the symbology of multiple layers at once.


texts: Altering linking arrow width in QGIS EasyCustomLabeling plugin

texts: Is it possible to alter the thickness of the line linking the label to the symbol, which is generated using the QGIS EasyCustomLabeling plugin?  I would like to make the line thinner.


texts:
I'm using the standalone GeoWebCache to serve tiles from a remote GeoServer. My problem is that the polygon label is added to each one of the tiles served, instead of only once in the polygon centroid.

I found a post which discusses the issue:
http://osgeo-org.1560.n6.nabble.com/polygon-label-repeated-for-each-tile-td4995203.html

The first reply mentioned a possible solution:

  "All in all, I suggest to use a tile rendering engine (GeoWebCache, MapProxy, TileCache) anyway, instead of requesting small image from GeoServer and have the tile rendering engine do the tile slicing afterwards. You will have send fewer requests to GeoServer (1 large image instead of multiple small images), so this speeds up the overall tile cache creation time."

Problem is that I couldn't find how to do that by referring to the GeoWebCache documentation, and the above mentioned post doesn't explain the way to implement that.

I also found a post with an answer that links to the GeoWebCache "Tiled" documentation, but my code allready uses all the necessary attributes and still the label shows up multiple times:

var Layer_1874 = new OpenLayers.Layer.WMS(
    'Grundkort',
    '/wms10.ashx'
    , {
        format: 'image/png',
        srs: 'EPSG:25832',
        layers: 'ballerupkommune_grundkort_bk',
        tiled: true,
        tilesOrigin: '698804,6173460'
    }
    , {
        displayInLayerSwitcher: true,
        isBaseLayer: true,
        transitionEffect: 'resize',
        displayOutsideMaxExtent: true,
        visibility: false
    }
);

Anyone has an idea?


texts:
I am looking for shapefiles with ward boundaries for the following Indian cities: Mumbai (electoral wards), Hyderabad, Ahmedabad, Surat and Kolkata.

Does anyone has them or know where I can get them from?


texts: Adding XY coordinate to map in ArcGIS (not in dataset)?

texts:
Is there a quick and easy way to add an XY coordinate in ArcGIS outside of creating a table?

A coworker just wants to see what that point is on a topo background.

I just want to see where a particular point is located without creating any data for it i.e. Just stored in a map document.


texts: identify sidewalks only at one side of the roadway

texts:
I got a task which I don't know how to save some labor work.

So I have two shapefiles in hand. One is roadway shapefile and the other is sidewalk shapefile. My task is to identify where there is only one-side sidewalk. The problem is when I created the sidewalk file, I didn't record the direction. Even if I did, I still can't come up with a way to automatically identify the roadway parts with sidewalks on one side.Is there anyway to make the identification easier? The only thoughts I have right now is to go through the sidewalk file and pick out the one-side sidewalks one by one.

Thanks for any of your thoughts.

Here is a screenshot.The red line represents sidewalks and the blue line represents the center line of roadway.


texts:
I have a problem with geocoding addresses (about one thousand). I have a file (.csv) with a large number of addresses: street, house number, ZIP code, city. I used the mmqgis plug, but after geocoding start QGIS freezes. Does anyone has similar problem? Should I use another tool?

QGIS version 2.14, operating system Windows 10,

There are no errors, only QGIS freezes. I did this before with another data (diffrent numbers) and everything was ok.

Now it works for few addresses and it freezes, or not at all.


texts:
I'm looking for a shapefile with this kind of information, only at considerably higher detail. In North America, ecoregions are classified by levels and grouped by vegetation type at much higher detail than on the linked map (which is of similar detail as the 14 WWF biomes). Using that classification, the data I'm looking for should be around level II or level III. E.g. it should differentiate between tropical rainforest, subtropical rainforest, and mangrove forests, and also it should with at least moderate accuracy show the tree line and the snow line in alpine regions.

(I found a WorldGrids raster file at somewhat high detail but it has a couple of problems, e.g. it shows actual vegetation rather than biomes, and I'm really looking for a shapefile anyway.)


texts: pgRouting algorithm for shortest path problem with resource constraints

texts:
do you know how to get pgRouting to approximate solutions to the shortest path problem with resource constraints (SPPRC), an NP-hard problem?

  SPPRC seeks a shortest (cheapest, fastest) path in a directed graph with > arbitrary arc lengths (travel times, costs) from an origin node to a
  destination node subject to one or more resource constraints.
  For example, one might seek a path of minimum length from s to t subject to the constraints that

  the total travel time must not exceed some upper bound and/or
  the total amount of some good that has to be picked up at the vertices along the path be less than or equal to some capacity limit and/or
  if two vertices i and j are visited on a path, then i must be visited before j

Quoted from the Boost Graph Library at http://www.boost.org/doc/libs/1_54_0/libs/graph/doc/r_c_shortest_paths.html.

I need such an algorithm to compute low-stress bicycle routes through an Open Street Map network from one origin to many destinations.
The 'low-stress' criterion introduces the constraints.

If no algorithm exists, then is there one that can at least compute all paths of a bounded length/cost from one origin to many destinations?
Then i can inspect each resulting path to find the ones that satisfy the constraints.
Probably not efficient, but it might be good enough for my purposes.


texts: Openlayers 3 - how to override handleUpEvent in PinchZoom interaction?

texts:
In Openlayers 3 I've been unable to get my attempts at overriding the ol.interation.PinchZoom's handleUpEvent to work.  I've tried a few different ways, but but below is a basic example of what I've tried.

This is an attempt to replace the standard PinchZoom's handleUpEvent with a custom one, and set this as the only interaction for the map.  However, the event function is never called (never logs anything).

function handleUpEvent(evt) {
    console.log("Up event handler");
    return true; //  Stop drag
}

map = new ol.Map({
    layers: [],
    target: 'map',
    controls: controls,
    interactions: [new ol.interaction.PinchZoom({handleEvent: handleUpEvent})],
    view: new ol.View({projection: projCode})
});


texts: Python script adding same data over and over when run from ArcToolbox?

texts:
I'm fairly new to python scripting.

I've written a script that runs fine from PythonWin or the ArcMap Python window, but is having issues when I run it from ArcToolbox.

The script adds a feature dataset to a geodatabases and then takes a group of fcs (all named S_FEP from multiple geodatabases) and merges them into a single fc in that feature dataset. It then moves on and creates another feature dataset and merges another group of S_FEP fcs.

When it is run from ArcToolbox the subsequent feature datasets get created, but the first set of S_FEP fcs are merged and put in the new feature datasets instead of a new group of fcs. This only happens when running it from ArcToolbox.

# Import Modules
print "Importing Modules: arcpy, os & fnmatch"
arcpy.AddMessage("Importing Modules")
import arcpy
import os
import fnmatch

# Set Top Level Variables
topDir = r"C:\Test\August"
#topDir = arcpy.GetParameterAsText(0)

MonthYear = "082016"
#MonthYear = arcpy.GetParameterAsText(1)

workingFldr = "FEP_XS_Working"

# Set Top Environment
arcpy.env.workspace = topDir
print "Environment set to:", topDir
arcpy.AddMessage("Environment Set")

arcpy.env.overwriteOutput = True
print "Overwrite Outputs = True"
arcpy.AddMessage("Overwrite Outputs = True")

# Create FEP_XS_Working Folder in topDir
arcpy.CreateFolder_management(topDir, workingFldr)
print "Created FEP_XS_Working Folder"
arcpy.AddMessage("Created FEP_XS_Working Folder")

# Create IDNR_FEP_XS_MonthYear.gdb in FEP_XS_Working folder
arcpy.CreateFileGDB_management(os.path.join(topDir,os.path.basename(workingFldr)), "IDNR_FEP_XS_" + MonthYear + ".gdb")
workDB = os.path.join(topDir, os.path.basename(workingFldr),"IDNR_FEP_XS_" + MonthYear + ".gdb")
print "Created Database:", workDB
arcpy.AddMessage("Created Database")

# List folders in topDir
for basinfolder in os.listdir(topDir):
    if basinfolder.startswith('0'):
        print "Found Basin:", basinfolder
        arcpy.AddMessage("Found Basin")

        # Create Feature Dataset
        cs = arcpy.SpatialReference('NAD 1983 UTM Zone 16N')
        arcpy.CreateFeatureDataset_management(workDB, "HUC10_"+basinfolder, cs)
        featureDS = os.path.join(workDB, "HUC10_"+os.path.basename(basinfolder))
        print "Created Feature Dataset:", featureDS
        arcpy.AddMessage("Created Feature Dataset")

        # Create List of S_FEP Layers and Create Variable for their Projections
        FEPMergeLst = []
        for dirpath, dirnames, filenames in arcpy.da.Walk(topDir):
            for filename in fnmatch.filter(filenames, 'S_FEP'):
                FEPMergeLst.append(os.path.join(dirpath,filename))

        # Merge all FEP layers into a single layer
        if FEPMergeLst:
            outDir = featureDS
            outNm = "Flood_Elevation_Pts_DRN_Water_"+os.path.basename(basinfolder)

            outLyr = os.path.join(outDir, os.path.basename(outNm))
            print "outLyr =", outLyr

            arcpy.Merge_management(FEPMergeLst, outLyr)
            print "Completed S_FEP Merge"
            arcpy.AddMessage("Completed S_FEP Merge")


texts: Thresholding distance raster and vectorizing in Google Earth Engine

texts:
I am struggling to convert a raster to a vector using Google Earth Engine.
I have produced an image using the ee.Image.distance() function where each pixel represents that pixel's distance in meters from the nearest green space (pixels with NDVI >= .5). The end goal is to produce a vector layer of the areas farthest from green space (specifically, regions over 1500m from areas with an NDVI >= .5) that can be used as the foundation of an intersection analysis (for example, what zip code areas overlap with these regions?).
I isolated the regions with a distance to green space greater than or equal to 1500m, and attempted to convert that raster image to a vector using code from GEE's 'Raster to Vector Conversion' guide. The raster image that should serve as the base of the vector conversion appears correct--the zones (less than 1500m, and greater than or equal to 1500) appear correct when visualized. However, when attempting to visualize the vector boundaries, they appear quite wrong. See the image below, where vector boundaries are in green and the raster areas they are meant to encompass are grey.

So, in short, how do I modify my code below so that I can subsequently convert the regions above that threshold to a vector FeatureCollection? The code is below, and a link to it can be found here.
// ==================== NDVI ====================

// Import the Landsat 8 TOA image collection.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA');

// Pull the least cloudy LANDSAT 8 image between May and August 2020 that has coverage
//   across the five boroughs
var image =
  l8.filterBounds(nyc_bb)
    .filterDate('2020-05-01', '2020-08-31')
    .sort('CLOUD_COVER')
    .filterBounds(roi)
    .first()
    .clip(nyc_bb)
    //.map(maskClouds)
;

// Compute the Normalized Difference Vegetation Index (NDVI) for that image.
var nir = image.select('B5');
var red = image.select('B4');
var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI').clip(nyc_bb);

// ==================== Distance to Green Area ====================

// Create zones using an expression, display.
var ndvi_zonesExpression = ndvi.expression(
    "(b('NDVI') >= .5) ? 1" +
      ": 0"
).clipToCollection(nyc_bb);

// Set Kernel
var kernel = ee.Kernel.euclidean({radius: 5000,
                                  units: 'meters'});

// Generate Distance image
var dist = ndvi_zonesExpression.distance({
  kernel: kernel
  }).rename('dist');
Map.addLayer(dist, {min: 0, max: 2600, scale: 30, palette: palettes.matplotlib.inferno[7]}, 'Distance from Green Space (m)');

// Take spaces furthest from green space (over 1500m) and convert to vector
var zones = dist.gt(1500);
zones = zones.updateMask(zones.neq(0));

var vectors = zones.addBands(dist).reduceToVectors({
  geometry: nyc_bb,
  crs: dist.projection(),
  scale: 30,
  geometryType: 'polygon',
  eightConnected: false,
  labelProperty: 'dist',
  reducer: ee.Reducer.mean()
});

// Display the thresholds.
Map.addLayer(zones, {min: 1, max: 3, palette: ['grey']}, 'thresholded and masked raster');
Map.addLayer(vectors, {color: '00FF00'}, 'vector boundaries');


texts: Is there anywhere that I can obtain global UV data in either raster or vector format?

texts:
I'm looking for some UV data (currently with no spatial or temporal resolutions in mind, but ideally split into UV-A and UV-B) for the whole earth. Can anyone recommend a source for data like this or point me in the right direction?

Any help much appreciated.


texts: Is QGIS compatible with PostgreSQL 8.4 Advanced enterprise DB with Oracle compatibility?

texts:
Sorry for such a long question.
We are trialling EnterpriseDB Postgres Plus Advanced Server 8.4 with Oracle compatibility and PostGIS 1.5.1.
I can create a connection to the database in qgis and the test connection is successful. However when I actually connect to the database to view the available tables i get the error below.
I understand that the utl_smtp is a schema related to the oracle compatibility portion of the install, so I assumed that opting for the value "Only look in the geometry_colums table" would have solved the problem. Unfortunaley not!
Has anyone else encountered this? Or more importantly solved this?

Thanks
Andrew


texts: Spearman correlation between two rasters in R

texts:
I have two rasters of the same area. The only difference is that they were created using two different methods. I would like to compare them and see how correlated they are pixel by pixel.
Any ideas how to do this in R or how to do a spearman correlation.
I wanted to do something similar to this (http://www.hakimabdi.com/20150101/an-r-function-to-test-pixelwise-correlation-between-two-time-series-of-raster-data/) with the result presented as a third raster.


texts: Add OpenLayers layer from QGIS's Python console?

texts:
On QGIS 2.4.0-Chugiak the Openlayers plugin is found under Web menu rather than the Plugins menu. How can you use the Python console to add a Openlayers Google Maps layer?

The following code give us the error OpenLayers plugin not loaded.

try:
    olplugin = qgis.utils.plugins['openlayers']
    ol_gphyslayertype = olplugin.olLayerTypeRegistry.getById(0)
    olplugin.addLayer(ol_gphyslayertype)
except KeyError:
    print 'OpenLayers plugin not loaded.'


texts: Adding a Geopackage layer as a hardwired input to an algorithm in the QGIS graphical modeler

texts:
Update- August 2020-  you can now load vector geopackage layers directly in 3.14 https://github.com/qgis/QGIS/issues/30827
I have a geopackage layer that I want to always be used in an algorithm in the QGIS 3.6 graphical modeler. When I use the triple dots to navigate to the geopackage it doesn't give me the option to choose a layer in the geopackage. Is there syntax I can use in the file path to identify the layer? In the example screenshot I want to use the layer AOI_Poly2 from the geopackage AOI as the mask layer. It seems like it will work but it chooses the first polygon layer in the gpkg (AOI_poly). I have tried appending  "|layername=AOI_poly2" to the file path but the dialog won't let me save. I am using windows 10.


texts: Which type of slash to use in path - QGIS 3.8 with Python

texts:
I'm a beginner using Python. I had a problem loading a shapefile in QGIS 3.8 using Python, while loading a rasterfile worked perfectly.

I found the solution here: Can't load a shapefile with Python Console in QGIS?

I loaded the files like this:

path = 'D:/Users/Me/data/'
rlayer = iface.addRasterLayer(path + "/landsuit.tif", "basemap")
vlayer = iface.addVectorLayer(path + "/ne_50m_urban_areas.shp", "urban", "ogr")

Loading the vlayer like this did not work, while loading the raster layer did. When I change the vlayer path to iface.addVectorLayer(path + "\\ne_50m_urban_areas.shp", "urban", "ogr"), thus replacing "/" with "\", the code works.

Can anyone explain to me why loading the vector layer with one single forwardslash does not work, while it does work for rasterfiles? In general: When is it better to use forwardslashes and in which cases should I choose double backwardslashes? Searching the Internet I have not found any satisfying solution.

I am using Windows 10 and QGIS 3.8.


texts: Installed Linux Ubuntu version 3.16.3 and want to split a state parcels file (95 counties), but when I click on the Vector menu item nothing happens. Have used this feature successfully before in an earlier version (don't remember which).  Also an earlier version had the ability to split out just one county by an option to specify which one, which seems to have now disappeared? Have reloaded v3.16 and problem still exists (can't seem to install 3.18)


texts: ArcGIS. Coordinates conversion from WGS84 to Stereo70 is not precise

texts:
So I want to convert coordinates from WGS84 to Stereo 70 and it is not working precise.
For example, I have the coordinates (20.2600, 43.6200) and in Stereo 70 should look like (448469.26,325966.67) but when I convert it, it gives me an error.

My line of code looks like:

Log.e("stereo70", " " + GeometryEngine.project(24.3512, 44.4316, SpatialReference.create(31700)));

Output: 448464.6330154536, 325965.9619150401
Which is about 30m away from the original point.

I have also tried a code provided in the answer:

        Point p = new Point(24.3512, 44.4316);
        Geometry projectedGeom = GeometryEngine.project(p, SpatialReference.create(SpatialReference.WKID_WGS84), SpatialReference.create(31700), GeographicTransformation.create(15995));
        Point projectedP = (Point) projectedGeom;
        Log.e("test", "X:" + projectedP.getX() + " Y:" + projectedP.getY());

Output: 448222.593914864 Y:325897.33784212614
Which is an error of about 300m.

How can I fix that?


texts: Showing filter options on embedded map of Carto?

texts:
With Cartodb, is it possible to show the filter options on the embedded map?

The explain, I have generated a visualisation and have shared, but I want people accessing the visualisation to be able to filter the results.


texts: Loops for calculating vegetation indices for multiple temporal images and each image have multiple bands

texts:
I have 64 TIFF files of different dates, each date has 4 bands (123 and 4 for B, G R and NIR). I am trying to develop a Vegetation index for each date image. This is how the images in my folder are

[1] "20171106.tif"     "20171106.tif.ovr" "20171113.tif"     "20171113.tif.ovr" "20171121.tif"     "20171121.tif.ovr" ........................................
 [64] "20171208.tif"

I used this command to process one single date image to get my index

st <- stack("20171106.tif")
br <- brick(st)
gcvi20171106 <- (br[[4]]/br[[2]])-1
plot(gcvi20171106)

But if I have to do it 64 times by manually changing the names its tedious.

I tried this loop

setwd("D:\\AP\\all_images\\")
lst <- list.files()
lngth <- length(lst)
rlist <- lapply(lst, raster)
for (i in 1:64){
  stk <- stack(rlist[i])
  brik[i] <- brick(stk[i])
  gcvi[i] <- (brik[[4]]/brik[[2]])-1
  }

I am new to loops. How can I do this?


texts:
I'm creating a video with the Time Manager Plugin in QGIS 3.11.

How do I draw data from an attribute for a tugboats engine power in a static location on the screen?

As it is it currently follows the boat as it moves along with the ship.  Once I start integrating more vessels it's going to get confusing doing it this way though.  Would prefer the data just display on the right hand side of the screen.


texts: Import raster to postgis and add a foreign key constraint

texts:
I use raster2pgsql to import raster files into a table with this command line:

C:\Program Files\PostgreSQL\9.3\bin>raster2pgsql -s 4326 -M -F -a c:\temp\mun.tif public.basqui_raster | psql -U postgres -d basqui -h localhost -p 5433

However I would like to add a foreign key constraint to each row during the import. Is it possible to do that when importing raster with raster2pgsql?

What could be the best way to do it?

Edit: I would prefer not use the filename (-F) in a second step to create the foreign key. I would prefer to do it during the import.


texts: Working examples of Leaflet/Polymaps with GeoJSON tiles generated from Tilestache using Shapefiles?

texts:
I'm looking to create choropleths using GeoJSON tiles generated from ESRI Shapefiles by TileStache(or similar), in either Leaflet or Polymaps. So far reading and drawing GeoJSON from both Leaflet and Polymaps is really easy, now what I'm looking for are working examples of Tilestache or anything like it generating GeoJSON from ESRI Shapefiles. Looking at the examples on the site, it seems to be pretty straightforward(just some config settings), but I'm wondering if there are any existing implementations out there that I could take a look at to see if I'm heading in the right direction. Google has found me nothing thus far.

Additionally, if anyone knows of a solution for generating Choropleths, or any other input regarding Choropleth generation, that would be fine as well. Really, any input here would be appreciated. Thanks!


texts: Convert a feature collection to image collection

texts:
I have created a feature collection of burned area in the Amazon according to 1-degree grid cells over 20 years. Here is a link to the collection - https://code.earthengine.google.com/?asset=users/stavo/FCburn

In this collection, features are the grid cells, columns names are dates.

In order to create a time-lapse video, I would like to create an image collection where each image is a representation of all features and 1 column(month).

I manage to create the wanted image for 1 month with the following code:

var Jan2015 = FCburn.select('2005_01_01_BurnDate')
                    .reduceToImage(['2005_01_01_BurnDate'], ee.Reducer.first())
                    .clip(borderB)

However, I can't manage to figure out how I would iterate through each column in the collection to create an image collection. I could not find an example of a function that iterates through columns rather than features in a feature collection.
How would you go about solving this issue?


texts: In a check-in/check-out syncronization setup, can a parent replica still be edited when child checked out?

texts:
In a check-in/check-out synchronization setup, can the parent replica geodatabase be edited while the child geodatabase is checked out and edited, and then when the child checks in their changes, have them reconciled against the parent changes?

I am trying to figure out what sync relation I can accomplish when limited to basic level license for the editing clients, but still allow the parent master to get changes while children are editing.  On that note, can a check-in/check-out have more than one child checked out at a time?

The main situation is that I am working with a small municipality who only have a Standard level license, but have three field crew tablets that need to be out in the field editing features simultaniously, the tablets connecting to an online version of the GDB at AmigoCloud.com which itself syncs back to the ESRI GDB.  But back in the office on the standard license I would like to continue working with the data while the field crews are out, and at end of day when they come back and sync I would like to sequentially sync each child with the parent and resolve any conflicts, but i wouldn't need to sync my changes to the tablets as my sync to the AmigoCloud child would migrate those edits.

So would the Standard license have to put out one child GDB, that child GDB is the one that syncs to AmigoCloud.com, and once all the field crews changes are into the AmigoCloud.com GDB that one is synced to the master and any conflicts are resolved, or with check-in/check-out is there no conflict resolution and one set of data is favoured over another?

I'd likely need to review the conflicts to settle them manually.  I guess after I am done syncing the AmigoCloud edits back to master GDB I would have to create a new child replica for AmigoCloud, so not sure this model would work given their sync.


texts: Interpolating values for NaN pixels using ENVI?

texts:
I'm currently using MODIS data and am running into problems with NaN values in the final product caused by clouds.  I'm interesting in a simple way to spatially interpolate (fill in) these NaN values.

Is there a menu function (in ENVI) that can easily/simply do this without delving into the intricacies IDL?


texts: Network service area mapping over open ground

texts:
Given a timestamped start location, a timestamped end location and a known average speed is there any way I can map everywhere someone (or something) could have possibly been between the two timestamps if they're able to travel in any absolutely direction in an area with no roads/paths etc?

The specific question here revolves around wildlife and camera traps. I have a dataset that includes details of when known individual animals were recorded by camera traps at known locations over a period of time. For example Elephant A was photographed at one location and then next seen five miles away fifteen hours later. Assuming an average speed of 3 MPH the elephant could have travelled from point A to B in a little more than an hour and a half. Being an elephant it's obviously wandered in a roundabout way (and probably stopped several times) in order to account for the extra 13+ hours between sightings.

What I'd like to do is to draw a shape covering the furthest extent in every direction that elephant could possibly have travelled in 15 hours and still made it to point B "on time". Obviously this will produce an exaggerated estimate since it's very unlikely any animal will have travelled constantly between sightings but I want to produce an overall idea of where they could have been which can (possibly) be narrowed down in future with more cameras.

Is this even vaguely possible?

I'm currently using QGIS and BigQuery's GIS tools but I'm willing to explore other solutions.


texts: How to randomly position coordinates within a containing feature in QGIS?

texts:
I have geocoded address data, and I am wishing to randomly position it within the census tracts that they are in to obfuscate... So the points are still within the feature but randomly move around so as to not be exact address locations, and to spread out some of the points that are all aggregated to a single point.

Is there a way to do this in QGIS?

One additional thing I need to point out, is that the geocoded data has additional data attached to it that I would like to preserve.


texts: Field Calculator Runs In IDLE (2.6.5) but not ArcGIS 10 SP3

texts:
just wondering if anyone else is having issues with the field calculator in ArcGIS 10 SP3. In the following code I do a pretty simple field calculation with a variable (an array) as one of the parameters in the expression.

exp1 = "Feature = 16"
rows = arcpy.SearchCursor("ObservationMaster", exp1)
li = []
for row in rows:
    ObsID = int(row.getValue("ObsID"))
    li.insert(len(li), ObsID)
if li != []:
    codeblock = """def chStat(stat, ID, li):
    if int(ID) in li:
        return 5
    else:
        return int(stat)"""

exp = """chStat(str(!ManagementStatus!), str(!ObsID!), li)"""

arcpy.CalculateField_management("ObservationDetail", "ManagementStatus", exp, "PYTHON", codeblock)

When I run it through IDLE (Python 2.6.5), there are no problems and it runs successfully. However, if I create a script in a toolbox in ArcCatalog/ArcMap using this same code, I get the following error:

<class 'arcgisscripting.ExecuteError'>: ERROR 000539: Error running expression: chStat(str(1), str(1), li) <type 'exceptions.NameError'>: name 'li' is not defined
Failed to execute (CalculateField).

I know this is because it is not treating "li" as a variable, but it does when I run it in IDLE. Has anyone had this same problem and does anyone know a solution to get this working. I really need this script to run via ArcMap.

Thanks in advance


texts: Calculating Ha of specific fields in QGIS?

texts:
I use QGIS Bonn 3.2.0
So my question is how can i calculate Ha for specific fields i marked on picture.
Is there any tool and if yes which and how to use it?


texts:
I would like to add two new columns in the table (E, N), fill in the coordinates of a polygon centroid using PostGIS SQL code and have the column upadated as soon as I save the edit of a drawing in QGIS.

Obs: MyGeometry is already in UTM


texts: Changing group icon in legend with Python in QGIS 2.6?

texts:
How to change the group icon in legend via Python in QGIS 2.6?

At version 1.x I did it this way:

self.legendTree = self.iface.mainWindow().findChild(QDockWidget,"Legend").findChild(QTreeWidget)
currGrp=self.legendTree.currentItem()
currGrp.setIcon(0,self.setMyIcon("folder.png")


texts: Efficiently changing field names added by Extract Multi Values to Points tool in ArcPy for large dataset?

texts:
How can I change the name of the composite raster of multiple bands as below- I think this may be the solution to efficiently address this problem but you are welcome for further solution?

I have a composite raster, Proximity_Rasters_NoHNSN.tif, which is produced by vertically stacking several other rasters.I need to extract all the band-wise values into a point layer in arcmap.

Now I am using Extract Multi Values to Points tool to extract band-wise values into a point layer which will be used for further analysis. But the problem is- the field added to the point layer after running this tool is not quite trackable since this tool adds fields like b1_Proximity_Rasters_NoHNSN.tif , b2_Proximity_Rasters_NoHNSN.tif...just prepending the b and band number to the raster name as you see.So I think running this tool after changing the band names will give me a meaningful and trackable result.

The field names added by Extract Multi Values to Points tool are as below:

How I am solving my problem now:

Now I am using arcpy to solve this problem by creating a new field and copying corresponding value into it and deleting the field after copying, but this way is too much time-consuming since I am dealing with large data.

Added code since @PolyGeo was guiding:

#Change field names
ffields = sorted([i.name  for i in arcpy.ListFields(eq_smpled_pnts) if re.match(r"b\d_", i.name)])
for old_fld in ffields:
    indx = int(re.findall(r'^b(\d{1,})', old_fld)[0])-1
    new_fld = raster_bands_names[indx]
    arcpy.AddField_management(in_table=eq_smpled_pnts, field_name=new_fld, field_type="DOUBLE", field_precision="12", field_scale="", field_length="25", field_alias="", field_is_nullable="NULLABLE", field_is_required="NON_REQUIRED", field_domain="")
    arcpy.CalculateField_management(in_table=eq_smpled_pnts, field=new_fld, expression= "!%s!"%new_fld, expression_type="PYTHON", code_block="")
    arcpy.DeleteField_management(eq_smpled_pnts, new_fld)

Update:
After trying @GISGe suggestion I found that all work correctly but the band1 does not change the name as intended!!See the screenshot of the script run as below:

Script I used to rename each band:

#Change composite band name
def band_name_changer(image_stack_path,names_ordered_list):
    raster_describe_object = arcpy.Describe(image_stack_path)
    if len(raster_describe_object.children) != len(names_ordered_list):
        raise Exception("Input raster names length and number of bands in the stacked image is not equal.......")
    else:
        for bnd in enumerate(raster_describe_object.children):
            indx_bnd = bnd[0]
            input_path = bnd[1].catalogpath
            output_path = os.path.join(os.path.split(input_path)[0], names_ordered_list[indx_bnd])
            arcpy.Rename_management(in_data=input_path, out_data=output_path, data_type="RasterBand")

But the output is ok except the BAND_1??


texts: Creating grid references with virtual fields

texts:
I wish to create a virtual field that will automatically update when new attributes are added; the equation is for working out the 10 figure grid reference of a point shapefile in SpatiaLite:

floor(((round($x)/100000)-floor(round($x)/100000))*100000) || floor(((round($y)/100000)-floor(round($y)/100000))*100000)

However, some of the rows appear with 'NULL' values - when I tried making a new field without virtual field selected they all came with a value but then this would not update when new points were added.

How can I fix this?

Also, this SpatiaLite layer will be shared across many projects, therefore will the virtual field be updated within each project?


texts:
I need to find a way to summarize an attribute in ArcGIS Online. I know it is a simple task in ArcGIS Pro where one just needs to look check the statistics and look at the properties section. However, I cannot find a way to do this analysis in ArcGIS Online.
How do I find the sum of the square footage for all entries while in an ArcGIS Online web app?


texts: I've got two tables, parcels and permits.  I'd like to display permits on a map without having to copy the geometry data from parcels.  Is there a way to create a join in Mapviewer?  Do I have to create a view?


texts: Leaflet - how to make geojson layer transparent (not visible on map),

texts:
I'm trying to make a geojson layer in leaflet transparent (not visible on the map) . I require the layer to perform a search based on road name using searchcontrol only. I already have the roads information attached as part of overlayMaps and use in layercontrol.

Below is my code:


texts: Understanding law or etiquette when applying changes to open-source software, specifically QGIS?

texts:
I am interested to know if there are any laws or etiquette surrounding the use of people's third party plugins or QGIS itself. What I am especially interested in is when a user accesses your code and amends it for their own needs, or if someone could take the source code of QGIS and rename it for commercial use.

I am not planning to do any of these things and even if I did I don't have the skills - this question arised from a debate in the pub


texts: replace data source for all layers that point to shapefile

texts:
I'm trying to do something quite simple but its not working. I want my python script to change all of the layers that point to the same datasource (a point shapefile) to another point shapefile. replaceDataSource seems to be the methods but when I run this:

for lyr in arcpy.mapping.ListLayers(mxd):
    #print lyr.name
    if lyr.dataSource == r"X:\filepath\folder\points.shp":
            print lyr.dataSource
            lyr.replaceDataSource(r"\\main.glb.corp.local\E-GB$\Home\AB\4\J022414\Documents\ArcGIS", 'SHAPEFILE_WORKSPACE', "ismage_template_P.shp")
            print "data replaced"
mxd.saveACopy(mxd)

i get "Runtime error type 'exceptions.ValueError'>: Layer: Unexpected error" when I run this. It prints the data source but doesnt print "data replaced" so it looks like the problem is on the lyr.replacedatasource line

Has anyone had this problem before, or should I be using a different method?


texts: Adding CSV layer to QGIS when X and Y fields are the same

texts:
I would like to add the CSV file, where both Longitude and Latitude head have been defined in one cell as per below:

But QGIS can't add them, as both X and Y are deemed as the same.
I've tried with WKT geometry (Option 2) but it's still wrong. The layer is inappropriate.

Is there any way to parse the following CSV data in QGIS?
All the numbers inside have been allocated to the specific coordinates defined in the top row and the first column.
I need them populated in QGIS.


texts:
I am trying to use arctoolbox and python so that the user can enter in parameters. In my script I have the following lines

_StationInterval = arcpy.GetParameterAsText(0)

I just want the user to be able to enter in a value for the station interval. However I get the error:

<type 'exceptions.NameError'>: name '_StationInterval' is not defined
Failed to execute (Final3Test).

Before I used arctoolbox, I had the station interval set to a value in the script, i.e.
_StationInterval = 10000

Am I missing something in my script??


texts: st_buffer() not preserving the shape of the original polygons

texts:
I am currently trying to create buffers of 5km around polygons of the region of interest. However, every time I increase the radius from 500 metres to 5000 metres, the shape of the buffer is not similar to the shape of the original polygons. How do I preserve the shape as I scale up the radius of the buffer?

Code:

trial <- st_read("C:\\Users\\vr235\\Downloads\\issuewithbufferusingsfpackage\\TownVillCity1.shp")
trial <- st_transform(trial, 32644)

Select 3 polygons from your shapefile to try

sel <- c(1,5,14)
geom <- st_geometry(trial[sel,])
buf <- st_buffer(geom, dist = 5000)
plot(buf, border = 'red')
plot(geom, add = TRUE)

EDIT The question has been modified based on one of the comments.


texts:
I have GML file with Application Schema whose data I need to import into SQL Server Spatial. I have been creted database structure with all thing that I need in database, and now I need to fill it with GML data.

How is this possible?

EDIT:

Does the SQLServer has some option to insert GML itself, or I need to use some other software like FME?


texts: Raster format with rotated cells - no resampling

texts:
Is there any way to rotate a raster w/o resampling? I am aware that rasters are always 90 deg perpendicular grids, but is there any file format or methodology that allows for a diagonal raster - one where the grid cells are actually angled akin to a vector?

I'm a new user, so I can't post images to show an example, but I want a raster to be rotated 18 degrees without the pixelation around the edges that is caused by nearest neighbor. We want original values, so bicubic and cubic convoultion won't do.

Thanks!


texts: SaveToLayerFile does not properly save symbology

texts:
I'm trying to save a symbology from a layer, but any way I try to do it, it just gives me a unsymbolised feature shapefile when I try to add the newly created .lyr file.

Anyone have any idea why this is?

I've tried using the ArcToolbox tool and it works fine, but the python snippet taken from that is a replic of what's in the script that's not working.

import arcpy
arcpy.env.overwriteOutput = True

head, wsName = os.path.split(r"D:\RS_Data\Workspace\Garden_CentralandNE")
mxd = os.path.join(r"D:\RS_Data\Workspace\Garden_CentralandNE",wsName + ".mxd")
outputFolder = os.path.join(r"D:\RS_Data\Workspace\Garden_CentralandNE", "gisdb")

mxd = arcpy.mapping.MapDocument(mxd)  # get MXD
df = arcpy.mapping.ListDataFrames(mxd)[0]  # get data frame
lyr = arcpy.mapping.ListLayers(mxd, "Service*", df)[0]

outLyr = os.path.join(outputFolder,"Symbology.lyr")

arcpy.SaveToLayerFile_management(lyr,outLyr,'#','CURRENT')  # Saves the file but without symbology


texts: Finding Sentinel tile for specific Long/Lat coordinate

texts:
I'm building a tool that downloads images from "Sentinel-2 on AWS"

I have a trouble translating Long/Lat coordinates to a specific Sentinel tile. I know that Sentinel 2 uses MGRS and I tried to simply convert Long/Lat to MGRS and with low precision I get the S2A tile ID. And it works, but not always.

For example, Long/Lat 34.665,31.625 resolves to "36R XA" in MGRS, but the Sentinel tile 36RXA doesn't exist.

What would be the right way of determining Sentinel tile using Long/Lat coordinates?


texts:
I have Geospatial application which is using GXP, GeoServer, GeoExt and PostgisSQL. We would like to enhance the application by adding a user friendly Admin system, similar to what GeoNode has. This sytem would let us create and manage users and create projects specific to users.

Can someone give us a head start on how we can do this? Or any place we can download it and configure it according to our requirements.


texts: Exporting multiple shapefiles using "select by attribute"

texts: I have about 12 different land cover shapefiles (all have common attribute data) that I need to select only pieces of it out.  After I run the select by location on all 12 shapefiles, is there a way (python, model builder or tool) to batch export all of the selected polygons from all shapefiles out at once instead of exporting them out one at a time?


texts: Removing Identical Polylines in multiple feature classes

texts: I have ~35 topo line files which show change in topography over time. 75% of the lines are identical since there was no change over the time period. I am trying to extract the portions of each file which have changed and delete the lines that are identical. Is there a way to batch process this?


texts: Identify same entities across time period

texts:
I have a few election district boundaries layers across several different elections, and I am trying to ascertain whether/how much the district boundaries have changed over time.

They are labeled with an ID but there is no guarantee that features with the same ID are in fact the same entity.

My rough idea is that I should compare the shapes of the boundaries across years, and if their overlap is sufficiently large (say 90%), then I can consider them the same entity.

I tried to google for an answer but being a novice I guess I just don't know the right keywords to search and couldn't find a way to achieve it.

Here is a screenshot of what they look like:


texts: Create Lines between points using points ID

texts:
I am trying to recreate in QGIS a few diferent sets of constellations. I already have a layer with all the stars I need, with ID and coordinates (downloaded from here astronomyData.zip), and I also have a table with a line for each constellation, and the pairs of stars that need to be connected by ID (link here)
What I need is for some way to have the software look at each pair of IDs, find the coordinate for each star with the same ID the stars layer, and connect each pair with a line.
I feel like this should be easy, but I can't figure out a way.


texts: Google Earth Engine getThumbURL lifetime?

texts:
How long can I rely on the URL produced by ee.Image.getThumbURL() to exist?

I have looked through the documentation, but I cannt find details about their longevity.


texts: Print FeatureLayer labels in Esri JavaScript API 3.5 and earlier

texts:
I have a need to implement labels on features in a FeatureLayer in version 3.5 if Esri's Javascript API. The labels are from a field returned by a REST feature service. I can't move to 3.7 for various reasons at this time. I have tried using a TextSymbol but my map features just turn to the color of the TextSymbol and no text appears. I may be approaching this in he wrong manner, though. Below is the code I'm attempting to use for labeling with the featureLayer object being my instance of the FeatureLayer I'm adding to the map. Is there a different or proper way to accomplish this task?

featureLayer.on("graphic-add", function (evt) {
var labelColor = new Color([255, 0, 0, 0.25]);
var myLabel = new TextSymbol(evt.graphic.attributes["My Field Name"]);

myLabel.setColor(labelColor);
myLabel.font.setSize("14pt");

evt.graphic.setSymbol(myLabel);

//console.log(evt);

});

Thanks for any help that can be provided!


texts: dimension 1 : coordinate intervals are not constant ERROR

texts:
I am calculating metrics with lidRpackage with the command grid_metrics applied to a catalog

opt_chunk_buffer(ctg) <- 0
opt_chunk_size(ctg) <- 0
opt_filter(ctg)       <- "-drop_z_above 50 -drop_z_below -1"
opt_output_files(ctg) <- paste(wd, "metrics/*_metrics", sep = "")

#opt_stop_early(ctg) <- FALSE #can be used if neccesary

#.options
opt <- list(raster_alignment = sqrt(250))

#METRICS CALCULATION
get_metrics = function(chunk)
{
  las = readLAS(chunk)
  if (is.empty(las)) return(NULL)

  if (area(las) <  sqrt(250)) {return(NULL)}                    #to avoid very small lasfiles that might have too few points
  else {
    metrics <- grid_metrics(las, ~stdmetrics_z.lasR(Z), sqrt(250), start=c(-0.2026,0))
    return(metrics)
  }
}
metrics <- catalog_sapply(ctg, get_metrics, .options = opt)

and some tiles give me an error which is:

metrics <- catalog_sapply(ctg, get_metrics, .options = opt)
#> |======                           | 7%
#> suggested tolerance minimum: 1
#> |======                           | 8%
#> An error occurred when processing the chunk 145. Try to load this chunk with:
#>   chunk <- readRDS("C:\Users\jaca\AppData\Loca\Temp\Rtmpc51D7l/chunk145.rds")
#>   las <- readLAS(chunk)
#> dimension 1 : coordinate intervals are not constant

It has happened to me before and I discovered that the tiles that retrieve the error normaly are very small and/or are not wall-to-wall together to the adjacent tiles. Because these tiles are not really relevant (they appear on edges of study area(see image)) the way I manage to "solve" my problem was to make a filter to only apply the function to those tiles bigger than certain area or to indicate opt_stop_early(ctg) <- FALSE but some times the first approach is not enough and the second one is not the best as if there are other errors they will be compiled.
Can someone explain me what is happening and if there is a way of solving it?


texts: How to display labels from many-to-many relationship in ArcGIS 10 using VB.NET?

texts:
I have two tables in a geodatabase:

parcel; and
person

1 or + parcel, can have 1 or +  person (relationship is many-to-many)

Using VB.NET I can display all the parcel labelled with the "id" field (see picture on left).

What I want is to display the parcel with labels for the Names of their persons (see picture on right):

I think that I have to use join between two tables

I can get the name with this code

 Sub addNom()
        ' Open the origin and destination classes.
        Dim featureWorkspace As IFeatureWorkspace = CType(Sp, IFeatureWorkspace)
        Dim originClass As IFeatureClass = featureWorkspace.OpenFeatureClass("parcelle")
        Dim destinationClass As IObjectClass = CType(featureWorkspace.OpenTable("personne"), IObjectClass)

        ' Create a memory relationship class factory and open a relationship class with it.
        Dim relClassFactoryType As Type = Type.GetTypeFromProgID("esriGeodatabase.MemoryRelationshipClassFactory")
        Dim memRelClassFactory As IMemoryRelationshipClassFactory = CType(Activator.CreateInstance(relClassFactoryType), IMemoryRelationshipClassFactory)
        Dim memRelClass As IRelationshipClass = memRelClassFactory.Open("pll_personne", originClass, "id", _
                                                destinationClass, "id", "Contains", "Is in", esriRelCardinality.esriRelCardinalityOneToMany)

        ' Find the counties related to the state with an ObjectID of 1.
        Dim originFeature As IFeature = originClass.GetFeature(6)
        Dim relatedSet As ISet = memRelClass.GetObjectsRelatedToObject(originFeature)
        relatedSet.Reset()

        ' Iterate through the set while finding the names of the counties.
        Dim nameIndex As Integer = destinationClass.FindField("nom_fr")
        Dim setObject As Object = relatedSet.Next()
        While Not setObject Is Nothing
            Dim row As IRow = CType(setObject, IRow)
            Console.WriteLine("OID: {0}, NAME: {1}", row.OID, row.Value(nameIndex))
            setObject = relatedSet.Next()
        End While
    End Sub

The Question is:

How  can I display the labels using this result?