texts: Exporting attachments and ignoring attributes that don't have attachments using ArcPy?

texts:
I'm using the following python script to export photo attachments:

from arcpy import da
import os

attachTable = arcpy.GetParameterAsText(0)
origTable = arcpy.GetParameterAsText(1)
nameField = arcpy.GetParameterAsText(2)
fileLocation = arcpy.GetParameterAsText(3)
with da.SearchCursor(attachTable,['DATA','ATT_NAME','REL_GLOBALID']) as attachCursor:
    for attRow in attachCursor:
        binaryRep = attRow[0]
        fileName = attRow[1]
        relID = attRow[2]

with da.SearchCursor(origTable,['GlobalID', nameField]) as originCursor:
    for origRow in originCursor:
        origID = origRow[0]
        origName = origRow[1]
        if origID == relID:
            break
        open(fileLocation + os.sep + origName + ".jpg", 'wb').write(binaryRep.tobytes())
        del origRow

I get this error after running it:

  Runtime error

  Traceback (most recent call last):   File "", line 7, in
   IOError: [Errno 2] No such file or directory:
  u'W:\work\dry_creek_exp\images\Hwy 20/Crowley Rd (turn left).jpg'

... I believe this file doesn't exist because I didn't have a photo attached to it. My question is, does anyone know of a way to get python to ignore these files that don't have attachments?? I hoping when I save this code as a script and run it on bigger projects that it will help the run time.


texts: Selecting surrounding states from a given state

texts:
I have a PostGIS spatial enabled table with the states in the US. How can I write a query selecting all touching states to a named state. Like this :

Select * from States Where geom touches (Select geom from States Where name = 'Oklahoma')


texts: Does switching satellites by GPS receiver cause errors?

texts:
Suppose some GPS receiver getting information ("connected") from 4 satellites, and from some reason disconnected from one of the satellites, and instead connecting to new one, so in total still there are 4 satellites as sources for information.

Is this possible? And if so, does some accuracy error (UERE) occur?


texts: Define compound classes by intersection of two raster images in Google Earth Engine

texts:
I am trying to create ecoregions by the intersections of soil classes and land cover classification.

Ideally I would like to have something like "forest_Entisol" and "forest_Oxisol" or as well "grass_Oxisol".

I have attempted to multiply the two classified images like this:

//import study area
var area = ee.FeatureCollection('users/jonathanreith/KK')
Map.centerObject(area,9);
//////////////////////////////SPATIAL JOIN//////////////////////////////////
var soil      = ee.Image("OpenLandMap/SOL/SOL_GRTGROUP_USDA-SOILTAX_C/v01");
var globc     = ee.Image("ESA/GLOBCOVER_L4_200901_200912_V2_3");
var ecounit = soil.add(globc);
var ecounit = ee.Image.cat([soil, globc]);
var ecounit = ecounit.select(
    ['grtgroup', 'landcover', 'qa'], // old names
    ['soil', 'landcover', 'trash']  // new names
);
var ecounit = ecounit.select('soil').multiply(ecounit.select('landcover')).rename('ecounit');
print(ecounit);

It kind of works, but I have no idea what is behind each value. Furthermore, classes are potentially not unique, because the product of the two values could have the same number as a original class.

Is there a better way in GEE to intersect two raster images?

I already saw Intersect of raster image on Google Earth Engine


texts: Store global reference to a service/object

texts:
In developing with ArcMap, is there any way to register a service that can be used during the life-cycle of a document- for example to do something like:

C# .NET

mxApplication.RegisterService(myWebService);

then later-

mxApplication.GetService<WebService>().SaySomething("Hi");

So kind of using the Service Locator pattern or similar.


texts: ArcGIS tool for indicating changes in land use

texts:
I have a couple of shapefiles in which there is some information about land use (polygons) such as water, forests, soils, buildings etc.
Each of the shapefiles represents a certain year of the same area of interest.

I want to indicate the area which was the biggest change over the years.

Can you suggest me a procedure to do it in ArcGIS 10 ?


texts: Finding Intersection Points in Stream Network using QGIS?

texts:
I have a stream network and I want to find the red points as shown below.

I have tried using Q-chainage and Line Intersection feature in QGIS but it gives more than expected number of points. The red points 513 and 556 are expected intersection points of stream segment and point 518 is expected start/end point, but the rest points are not required.These extra points may occur because the poly line representing each stream segments is broken into multiple line segments represented by intermediate points which is expected to be single geometry in ideal case.

So, I have also tried merging the poly-lines based on "value" attribute (in PosgreSQL by ST_Union ) as it had common value for most of the stream segments. But then in some cases still the split of poly-line is not for each stream segment but spans out to other stream segments as well.

Please also point out if there are errors in pre-processing of data. I have generated the stream network from DEM with and without eliminating the depressions using r.fill.dir in QGIS . Is there any other processing required?
If not, please point out a method to find only start,end and intersection points in the given network.


texts:
Would it be possible to calculate a density metric in height intervals as in FUSION or Lastools?
e.g number of all returns between 0.2 and 3 meters/number of total returns in 1 meters cell.
I'm trying to compute using the function of the package stdmetrics_z(z) but I had an error with fast_table

myMetrics2 = function (z) {
  n = length(z)
  zmax = max(z)
  zmean = mean(z)
  if (zmax <= 0) {
    d = rep(0, 9)
  }
  else {
    breaks = seq(0.2, 25, 2.8)
    d = findInterval(z, breaks)
    d = fast_table(d, 10)
    d = d/sum(d) * 100
    d = cumsum(d)[1:9]
    d = as.list(d)
  }
  names(d) = paste0("zpcum", 1:9)
  metrics = c( d)
  return(metrics)
}

metric2=grid_metrics(ctgn.i, ~myMetrics2(z, res=1)
#> Processing [=========================================] 100% (1/1) eta:  0s
#> Error: could not find function "fast_table"


texts: Converting polygon values to a vector grid using PostGIS

texts:
What is the best way of converting a value from irregularly shaped poylgons to a polygon grid using PostGIS? I have a set of polygons that look like this:

and I would like to overlay a grid on these polygons, and assign values from the polygon to the grid value, proportional to the area of the polygon that intersects the grid. For example, where a grid intersects more than one polygon, I want to use the fraction of each polygon area that intersects with the gridcell to proportionally calculate the value.

My approach for doing this is as follows. Assuming each polygon has a value X_i and an area A_i, I calculate the density using D_i = X_i / A_i. Then, for each gridcell, I multiply the fraction of the polygon area that intersects that gridcell by the appropriate density value. I think this approach is satisfactory but I have not figured out a way of abstracting it so that it is easier for me to do this regularly.

I am curious to learn if there is a better way to do this, if there is a term for this process and if there is a function that exists to perform a calculation of this type.


texts: EPSG:4326, EPSG:3857, JOSM and unprojections

texts:
It is said here that EPSG:4326 is not a projection.
In my mind, a projection is a transformation of the surface of the Earth producing a flat image, notably a bitmap.
We know that WGS84 (World Geodesic System '84) is a model to reference the Earth surface.
It's an angular coordinate system and not a projection.

When I request an image from a WMS server and I use EPSG:4326, I receive a bitmap image. Hence, in my mind, EPSG:4326 is both a projection and a coordinate system used to locate the image (its dots).  The projection is Mercator Equidistant Cylindrical  and the coordinates are WSG84.

JOSM is said to be unable to reproject a bitmap, that is, to transform a bitmap image from one projection to another. However, if I choose EPSG:3857 as the projection JOSM uses to display the map "at the screen" (in a window), JOSM is quite able to display 4326 from the server (but not other projections). Hence, the image must be "the same", more exactly "not to be reprojected".
This seems incoherent with what you have just answered (4326 and 3857 displaying differently). So what's going on there?

Have I got that right and is EPSG:4326 a projection or not? (in bold on iant's request)
BTW, I read this, seeming to say that an EPSG is both projection and coordinates:
"Les codes EPSG constituent une liste des systèmes de coordonnées géoéréférencées de projection"


texts:
In server status page I can see JAI Memory usage is always 0KB.

I already enabled Tile Recycling.
Is it related to the JAI-EXT? I activated it using -Dorg.geotools.coverage.jaiext.enabled=true in startup script.

Why is it always 0KB? How much is the server performance impact and what can I do to force/make my server use JAI memory?


texts: Creating coverage file (vector) from raster data and indexing it to original files?

texts:
I have a large number of LiDAR files (in .asc format) that I'm wanting to catalogue so it is easier for people to find the files they want.

How can I take LiDAR files I have, and create a .shp file (ideally polygons) of the areas covered.

My end goal is to have a .shp file people can open on top of a map, check whether or not their location has been covered by LiDAR and then identify the file name that they're after.


texts: stars package in R. How to add time dimension in a stars object?

texts:
I have a raster stack (my_stack) and I converted it into a stars object:

mystar <- st_as_stars(my_stack)

That results in:

stars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
    layer.1
 Min.   :  0.0000
 1st Qu.:  0.0000
 Median :  0.0000
 Mean   :  0.7897
 3rd Qu.:  0.0000
 Max.   :686.6100
dimension(s):
     from  to  offset delta                       refsys point               values
x       1 904  104010  1000 +proj=tmerc +lat_0=0 +lon...    NA                 NULL [x]
y       1 773 4623806 -1000 +proj=tmerc +lat_0=0 +lon...    NA                 NULL [y]
band    1  67      NA    NA                           NA    NA layer.1,...,layer.67

I also have a vector with Date objects.
How could I add this vector as time dimension in mystar object?


texts: Get weird result when calculating distance between two points, using latitude longitude?

texts:
I use algorithm (see: calculating distance between two points) to calculate distance between two points. but it get wrong result, see Java code:

/**
 * group1：set latitude to 89
 */
// set longitude1 = 0, longitude2 = 1
System.out.println(distance(89, 89, 0, 1, 0 , 0)); // 1941.0            result1

// set longitude1 = 0, longitude2 = 2
System.out.println(distance(89, 89, 0, 2, 0 , 0)); // 3881.0            result2

// set longitude1 = 0, longitude2 = 179
System.out.println(distance(89, 89, 0, 179, 0 , 0)); // 222381.0        result3

// set longitude1 = 0, longitude2 = 180
System.out.println(distance(89, 89, 0, 180, 0 , 0)); // 222390.0        result4

/**
 * group2：set latitude to 0
 */
// set longitude1 = 0, longitude2 = 1
System.out.println(distance(0, 0, 0, 1, 0 , 0)); // 111195.0            result5

// set longitude1 = 0, longitude2 = 2
System.out.println(distance(0, 0, 0, 2, 0 , 0)); // 222390.0            result6

// set longitude1 = 0, longitude2 = 179
System.out.println(distance(0, 0, 0, 179, 0 , 0)); // 19903892          result7

// set longitude1 = 0, longitude2 = 180
System.out.println(distance(0, 0, 0, 180, 0 , 0)); // 20015087          result8

In group1:

a = result2 - result1 = 1940

b = result4 - result3 = 9

In group2:

c = result6 - result5 = 111195

d = result8 - result7 = 111195

It's very strange. c = d, group2 is right, but a ≠ b, group1 is wrong.

How do I solve this issue?

Is there other higher accuracy algorithm for calculate distance between two coordinate points?


texts:
I'm working with London OSM data (roads) with PostGIS and pgRouting.

After download the pbf file from Geofabrik and create the topology with osm2po, I've realized that some highways in London are missing (the M25 Highway around Greater London)
I've tried both pbf files: Britain Islands and Great Britain, and in both is missing, but if I search for this highway via Overpass Turbo, I can download it.

¿Has anybody any idea why is this road missing in this dataset? I thought that geofabrik datasets were daily extract from OSM Database.

As a solution, I've thought in manually create the road, and then calculate the topology. ¿Could this be a solution, or would I have to do something more in order to add it to all network?


texts: Import qgis.core in PyCharm: "No Module Named Pyqt5.Qtcore"

texts:
I would like to work on Python script with tools from qgis in a standalone software (like Pycharm or IDLE). I am quite new with Python programming and with setting up a work environnement.

In order to create the link between Python and the qgis module, i made the Batch file like this:

@echo off
SET OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
call "%OSGEO4W_ROOT%"\apps\grass\grass76\etc\env.bat
@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
path %PATH%;%OSGEO4W_ROOT%\apps\grass\grass76\lib
path %PATH%;C:\OSGeo4W64\apps\Qt5\bin
path %PATH%;C:\OSGeo4W64\apps\Python37\Scripts

set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python
set PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37

set PATH=C:\Program Files\Git\bin;%PATH%
start "PyCharm aware of Quantum GIS" /B "C:\Program Files\JetBrains\PyCharm\bin\pycharm64.exe"
cmd.exe

After I run this file, I obtain a cmd.exe window in which I try to import qgis.core and PyQt5 to see if it works correctly:

It looks like it's working. Perfect! Then, this .batch file opens PyCharm in which I try to import qgis.core in a new script, after stting up that the Python Interpreter is the on from QGIS (direction :C:\OSGeo4W64\apps\Python37\python.exe)

But when I try to import qgis.core, I get this message:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Program Files\JetBrains\PyCharm\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\OSGEO4~1\apps\qgis\python\qgis\__init__.py", line 70, in <module>
    from qgis.PyQt import QtCore
  File "C:\Program Files\JetBrains\PyCharm\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\OSGEO4~1\apps\qgis\python\qgis\PyQt\QtCore.py", line 24, in <module>
    from PyQt5.QtCore import *
  File "C:\Program Files\JetBrains\PyCharm\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
ModuleNotFoundError: No module named 'PyQt5.QtCore'

However, when I try "import PyQt5", it looks like it's working, but still importing qgis.core doesn't work and leave me with the same message.

Does anyone have an idea of how I can fix this?

I just want to be able to script tools with the qgis.core module but I'm out of ideas.


texts:
I have about 7000 records in a feature layer. It was all fine but when I updated the features of record, suddenly on ArcGIS Online when I open the "Data" tab first it shows data loading (7000 records) and suddenly after that it shows "data error" and I don't see any previous data.
What can I do?


texts: How to check if two qgsPoints are equal in PyQgis

texts:
I have two points:

a=QgsGeometry.fromPoint(QgsPoint(227739,908343))
b=QgsGeometry.fromPoint(QgsPoint(227739,908343))

but when I run:

a==b

it returns False.
What am I doing wrong?

Update:

it returns True when I only use:

a=(QgsPoint(227739,908343))
b=(QgsPoint(227739,908343))


texts: Could not create tool plugin with ptype: gxp_zoomtolay

texts:
I am new in OpenGeo Suite Library and SDK client. I am following the tutorial for the tools:

http://suite.opengeo.org/opengeo-docs/webapps/sdkintro/viewer/featureinfo.html
http://suite.opengeo.org/opengeo-docs/webapps/sdkintro/viewer/zoomtolayerextent.html

But I get this error when I try to refresh my application:

  Error: Could not create tool plugin with ptype: gxp_wmsgetfeatureinfo
  src/ext/ext-all.js Line 7

My code is not long as I just started the tutorial and I havent modified much.
What I am doing wrong. Thanks in advance for your guidance.

/**
 * Add all your dependencies here.
 *
 * @require widgets/Viewer.js
 * @require plugins/LayerTree.js
 * @require plugins/OLSource.js
 * @require plugins/OSMSource.js
 * @require plugins/WMSCSource.js
 * @require plugins/ZoomToExtent.js
 * @require plugins/NavigationHistory.js
 * @require plugins/Zoom.js
 * @require plugins/AddLayers.js
 * @require plugins/RemoveLayer.js
 * @require RowExpander.js
 * @require overrides/override-ext-ajax.js
 * @require plugins/ZoomToLayerExtent.js
 * @require plugins/WMSGetFeatureInfo.js
 **/

var app = new gxp.Viewer({
    proxy: "<PROXY_URL>",
    portalConfig: {
        layout: "border",
        region: "center",
        // by configuring items here, we don't need to configure portalItems
        // and save a wrapping container
        items: [{
            id: "centerpanel",
            xtype: "panel",
            layout: "fit",
            region: "center",
            border: false,
            items: ["mymap"]
        }, {
            id: "westpanel",
            xtype: "container",
            layout: "fit",
            region: "west",
            width: 200
        }],
        bbar: {id: "mybbar"}
    },

    // configuration of all tool plugins for this application
    tools: [{
        ptype: "gxp_layertree",
        outputConfig: {
            id: "tree",
            border: true,
            tbar: [] // we will add buttons to "tree.bbar" later
        },
        outputTarget: "westpanel"
    }, {
        ptype: "gxp_addlayers",
        actionTarget: "tree.tbar"
    }, {
        ptype: "gxp_removelayer",
        actionTarget: ["tree.tbar", "tree.contextMenu"]
    }, {
        ptype: "gxp_zoomtoextent",
        actionTarget: "map.tbar"
    }, {
        ptype: "gxp_zoom",
        actionTarget: "map.tbar"
    }, {
        ptype: "gxp_navigationhistory",
        actionTarget: "map.tbar"
    },{
        ptype: "gxp_wmsgetfeatureinfo",
        outputConfig: {
        width: 400
        },
        actionTarget: {
            target: "map.tbar",
            index: 1
        }
     },{
        ptype: "gxp_zoomtolayerextent",
        actionTarget: ["tree.tbar", "tree.contextMenu"]
    }],

    // layer sources
    sources: {
        local: {
            ptype: "gxp_wmscsource",
            url: "/geoserver/wms",
            version: "1.1.1"
        },
        osm: {
        ptype: "gxp_osmsource"
        }
    },

    // map and layers
    map: {
        id: "mymap", // id needed to reference map in portalConfig above
        title: "Map",
        projection: "EPSG:900913",
        center: [-10764594.758211, 4523072.3184791],
        zoom: 3,
        layers: [{
            source: "osm",
            name: "mapnik",
            group: "background"
        }, {
            source: "local",
            name: "usa:states",
            selected: true
        }],
        items: [{
            xtype: "gx_zoomslider",
            vertical: true,
            height: 100
        }]
    }

 });


texts: grass algorithm changes integer to fixed length

texts:
using grass v.clean in the graphical modeller it changes all integer values from the original file to int with length 9

Is this a standard behaviour?  Can this be avoided?

The details of the grass algorithm are the following


texts: Combine slope info with Landsat for classification in Google Earth Engine

texts:
I wanted to add slope data as a band to filtered Landsat data.
I used ee.Image.cat to "merge" two images. I can see added band in a merged image, but when I perform the classification (which worked before adding slope data), I got error: Property 'sr_atmos_opacity' of feature '1_0_0' is missing.
Probably, I have to find another way how to add slope band to my data, but I can't figure out how to do that.
This is the classification that work (only Landsat, without slope): https://code.earthengine.google.com/bac39dc10202d7ca2e4ce37ed6b458c5
This script doesn't work (with slope):
https://code.earthengine.google.com/9b216bd7b78925e44e306b9e051aa73b


texts: Counting the number of layers in a GeoPackage

texts:
I'm trying to figure out a best practice for determining what are the names of layers in a GeoPackage and counting the total number of them.

I've been able to list all the layers of a GeoPackage using gdal info mypkg.gpkg but it only lists the layers, not count them.

The solution I need is not necessarily for scripting but simply for verifying whether my processes have successfully run. If there is a method to check via bash, QGIS3, Python, that would be very helpful.

EDIT: I should specify that the contents would include not only vector layers but raster layers as well.


texts:
I managed to learn the fundamentals of Django, GeoDjango and know OpenLayers for a while. Now I realize a project with a lot of objects with point geometries. To make the user map fast and responsive, I like to transfer only the points that are in the bounding box of the current view.

OpenLayers has on part of the solution on clientside solution with the BBOX strategy (new OpenLayers.Strategy.BBOX()) which appends the extends to the data URL e.g. /api/node?bbox=13.898900146484,51.04243240501,13.341099853516,51.117559096132

Now my question is, how I can deal with that on server side, without needing to reimplement attribute parsing and filtering on my own?
If there is no ready solution, might anybody please give me a rough snippet on what to do on django side?


texts: Performing find and replace in QGIS field calculator?

texts:
I'm trying to replace selected string values in an attribute table, eg. "Text1" for "Text2" (QGIS 1.8.0), patriculary for those rows with NULL values. I have been trying to use the replace expression (and any other that I have found) with no luck.

Are there any wildcard characters that I may be able to use in future?


texts: This may be trivial question. I am new to GIS. Could any of you can brief me the difference between the feature and Geometry.


texts: Conditionally require question in Survey123 if an image is submitted

texts: In Survey123 Connect (Desktop), is it possible to conditionally require a question only if the field worker has submitted a photo of an asset?


texts: add attributes to newly created shapefile from existing shapefile using arcpy

texts:
I have created a new shapefile and want to add shape attribute from another existing shapefile. It will be a recursive process with update of shape attributes of new shapefile (ofc) more than 2 times. Any help is appreciative. Following is what I am trying to make work:

Import arcpy
ofc = "C:/outFC.shp"
ifc = "C:/inFC.shp"
fields = arcpy.ListFields(ifc,"*","ALL")
for f in fields:
    shape = f.Shape
    query = '"Shape" = \'' + str(shape) + '\''
    cur = arcpy.UpdateCursor(ofc,query)

del f, fields

Diagram: I want to get FID (0) from first table and add it to second table of new shapefile and only need FID, Shape, and id fields in second table.


texts: Applying SetNull to multiple values at same time (array) gives Syntax error?

texts:
This is not working

I've got this error:

I've tried several combinations:

SetNull("rocksoilver12"==10 OR "rocksoilver12"==14, "rocksoilver12")

SetNull("rocksoilver12"==10 or "rocksoilver12"==14, "rocksoilver12")

SetNull("rocksoilver12"==10 OR "rocksoilver12"==14, "rocksoilver12")

SetNull("rocksoilver12" IN (10,14), "rocksoilver12")

SetNull("rocksoilver12" IN [10,14], "rocksoilver12")

SetNull("rocksoilver12" in [10,14], "rocksoilver12")

This is the error message that I've got:

None of them worked :/

Now, this works fine:

SetNull("rocksoilver12"==10, "rocksoilver12")

How can I use SetNull to delete multiple cell values at a time?


texts: Seeking alternative for ERDAS IMAGINE AutoSync?

texts:
I had the possibility to work with a Demo-Version of the ERDAS IMAGINE module AutoSync.
I have different geo-rasters (i.e. geoTIFFs) produced by UAV photography and which have to be matched the best way possible in order to carry out some analysis. If possible, even inner distortions should be corrected in order to match the reference raster.

The AutoSync module above does the trick, but as a small startup company, we are not able to pay the prices the company asks for the license.

Consequently, I am searching for an alternative to this software.

Do you know a good alternative?

Ideally it would be open source, but even a proprietary software would be OK, if the price is accessible.


texts: How to interpolate while providing clipping geometry? (ArcGIS 10)

texts: Is there a way to interpolate (IDW) while providing clipping geometry? My objective is that I want to interpolate temperature data but I do not want it to interpolate the entire extent of the point file but rather fit the output to the extent of a particular polygon. Thanks.


texts: Get Corresponding OSM Node IDS from list of Lat,Long (CSV): Map Matching

texts: I have a list of (lat,lng) coordinates. I want to map match them and get their corresponding OSM node ids? How can i do that using OSRM or any other map matching algorithm?


texts:
I am a novice, self-taught QGIS user (using v3.10.0) and I would like to integrate the National Library of Scotland (NLS) Tileserver API into QGIS, but I just can't crack it.

Tile server is https://nls.tileserver.com/

I am using CRS EPSG:27700 - OSGB 1936 / British National Grid - Projected

I would like to be able to switch to the API maps at varying zoom levels.


texts:
I would like to use the GEOS R-tree implementation to store the items of a growing structure.
However, I also need to query from the tree during the growth process to make sure that growing substructures do not overlap.

This means that I need to alternatively query from and insert to the tree.
However, qerying from GEOS R-tree triggers the build function, and inserting to a "built" tree is forbidden.

My question is therefore "is there a way around this?".
From my (limited) understanding of R-trees, this should be possible: there are supposed to support this kind of dynamic update without having to be rebuilt from scratch, so how can I do that with GEOS?
I've seen nothing in the C-API, and not much more than the cause of my problem:

// Cannot insert items into an STR packed R-tree after it has been built
assert(!built);


texts: ArcGIS Raster Calculator modifying bit depth of outputs?

texts:
I am trying to multiply rasters (converted from polygons and polylines) whose values I have reclassified into one class respectively with a value of 1. Then I input the weighted average equation in Raster Calculator, however I end up repeatedly with the following, devoid of statistics or tabular values and without anything showing up on screen:

I'm aware that the numbers I'm observing are indicative of the full value range for 32 bit depth (-2147483648 to 2147483647), although my DEM raster specification are 3186, 2073 (columns, rows) with a resolution of 65.15 X 65.15 meters, and 16 bits (and yes, I've made sure that my model environment's processing extent is specified to my DEM raster's specifications; and so are all the layers I'm processing).

Alternatively I tried using weighted overlay, which brings me the same issue but with a different range:

The value corresponds to the full value 8-bit range.

So clearly something is happening to the bit values when I'm trying to process my rasters together. I thought initially it was an issue with my output values being too large, but I doubt that's the case. What am I missing? Is it something to do with bit-depth?


texts: Mapinfo incorrect polyline (boundary) map to region conversion

texts:
How to solve issue regarding incorrect boundary to region conversion? I have boundary map plotted on cosmetic layer when I convert it to region it results in line with overlap to other region.
The two images are :

How do I resolve this issue?


texts: Applying scalars from a timeseries to raster classes in GRASS GIS (or python)

texts:
I have some yearly crop data layer rasters to which I would like to apply monthly crop-specific evapotranspiration depths from a text file. I could go month by month and do some raster math, but that seems like the clunky way to go about it. There are also several weather stations with different ET values for each crop, so my clunky expression would be something like

r.mapcalc expression = "ETdepths2015 =
    if(  CDL_2015@IDlandcover == 1 & stationfootprint@IDlandcover == 12, 0.3) +
    if(  CDL_2015@IDlandcover == 4 & stationfootprint@IDlandcover == 12, 0.4) +
    ..."

Problem is there are 50 unique crop class, 12 months in the year, 10 years, and 9 stations, so that would be about 54000 lines of mapcalc

I'm new to scripting in GRASS and I would like to use numpy to do the calculations without having to export the rasters, but I haven't gotten to far in that regard. I'll keep working on it, though if anyone has any hints, I'd be happy to hear them.

Update and some ideas

I was hoping use grass.script.array.array() to convert the rasters to numpy arrays and match crop codes, though I couldn't get that to work.

Here is a sketch of how I think I could still use r.mapcalc by building an expression iteratively.

# Loop ETI  years
for root,dirs,fnames in os.walk(rastpath):
    for fname in fnmatch.filter(fnames,'ETI_Idaho*'):
        ETIrast = fname  # String of raster map
        # TODO get year from raster map
        year = 2015
        # loop weather stations
        for site in d.keys():
            # loop through months and make a raster for each month
            for month in range(1,13):
                # Boolean array for grabbing specific row
                idx = dates == dt.date(year,month,1)
                # Loop through
                for crop in hdr[1:]:
                    # Get ET for a crop at site during month during year
                    ET = d[site][crop][idx]
                    # Mapcalc expression to create monthly raster
                    expr = '"ETdepths_%s_%s" = if( ETI_TV_voronoi@IDlandcover == %s & ETIrast == %s, %s ) + ' % (year, month,site,crop, ET)

And....I'm stuck on how to do that iterative string building!


texts: Getting a "square" of lat/lon for database query

texts:
So I want to get the square formed by a single point of lat/lons, my plan was to take that point's latitude, divide by 360, and multiply by the circumference of earth. Giving us the distance from the equator.

Now I will take 1/4 the circumference of earth, and subtract by the distance from the equator, giving us the height of triangle T. I can the find the base because I know the initial height, truncated height, and the original base.

Now using the base and height I can use trig to find the minimum and maximum longitude, and the latitude is easy because lets are evenly spaced. Does this sound right?

EDIT:

Let me simplify this question.

How can I find the distance between two points of longitude, at a single point of latitude?


texts: ATLAS automation in QGIS showing data from multiple layers (based on attribute from one of the layers)

texts:
I have a QGIS project with four layers: A, B, C and D.

I need to prepare an ATLAS automation based on layer A. On maps I want to see data from layers B, C and D (from A not necessarily) but the data needs to match to the attribute from table A selected for automation (the Coverage layer = Layer A, Filter with = values from 'Feature ID' column, Sort by = values from 'Feature ID')

Example:
Layer A - ATLAS automation based on values in attribute table from column: 'Feature ID' = 100, 101, 102 and so on.

Data from layers B, C, D visible on maps: features from layer B (C, D) where values in each attribute table from column 'Layer A Reference No' = 100, 101, 102 and so on.

I don't want to see the other features from layers B, C and D (where values are different than the ones selected for ATLAS automation (for example 87, 88 or 94)) for layer A - only the ones which match.

In the end I would see on the first map the data from layers B, C and D where attributes in the column 'Layer A Reference No' = 100 as this is the first map selected by ATLAS (based on Coverage layer A).

On the second map (where ATLAS selects value 101 from layer A) I would see only the data from layers B, C and D where attributes in the column 'Layer A Reference No' = 101.

Is this possible without Python?

At the moment the only thing I can do is to query layers B, C and D in order to see on the map the data where the value in the column 'Layer A Reference No' = 100 (for the first map), 101 (for the second map) etc. and then use layer A for setting location.

But I have to prepare 200 maps!


texts: PostGIS Zonal statistic (St_SummaryStats()) NULL problem

texts:
I try to assign each Polygon in my data set the maximum value that can be found within their areas. I have a raster and a collection of Polygons with each have a unique id (gid).

I uploaded the raster with:

raster2pgsql -d -I -C -M -F -t auto -s 32637 <Mypath> <Rastername> | psql -d <DBname> -U <Username>

When I run following code I get some results which appear to be reasonable for the segments, but I also get a lot of NULL values. When I run zonal statistics on the same geometry in QGIS (Raster Analysis > Zonal statistics) I get no NULL values while running GRASS v.rast.stats gives me exact the same output like in PostGIS, additional with the warning:

WARNING: Not all vector categories converted to raster. Converted 2128 of 3361.

Which correspondents also with the numbers of missing values in Postgis. I'm breaking my head whats causing this problem. I guess its related to the geometry created with ST_Buffer, wich somehow is not compatible with St_SummaryStats(). I also thought the problem is related to the tiling but even without tiling it occurs. ST_IsValid() returns for every element of my geometry TRUE, also is every element according to GeometryType() a Polygon.

  WITH blocks AS (
    SELECT row_number() OVER() AS gid,
           ST_Buffer(ST_Line_Substring(the_geom, 8*n/length,
             CASE
               WHEN 8*(n+1) < length THEN 8*(n+1)/length
               ELSE 1
             END), 4, 'endcap=flat join=mitre') AS the_geom
    FROM (SELECT ST_LineMerge(geom) AS the_geom,
              ST_Length(geom) As length
              FROM <LineDataset>
         ) AS t
    CROSS JOIN generate_series(0,CEIL(length/8)::int) AS n
    WHERE n*8/length < 1
  )
  SELECT gid,
         (ST_SummaryStats(St_Union(ST_Clip(rast, 1, the_geom, true)))).max AS a
  FROM <Rastername>, blocks
  WHERE ST_Intersects(rast, the_geom)
  GROUP BY gid;

The following picture shows the attribute table of the output of QGIS zonal statistics (left) and the code from above (right).

Some example data to recreate the error:

Some random Landsat data
Upload raster with: raster2pgsql -d -I -C -M -F -t 100x100 -s 3763 <somepath>\\Landsat8_L1TP_RGBN.tif test1 | psql -d <DBname> -U <Username>
Create sample line e.g.:

CREATE TABLE line AS
SELECT
ST_SetSRID(ST_MakeLine(ST_MakePoint(-28718,172580),    ST_MakePoint(-14012,182541)), 3763)::GEOMETRY(LINESTRING, 3763) AS geom

Run code from above after setting <Rastername> to test1 and <LineDataset> to line.


texts: How to select all tags in "Export OpenStreetMap topology to SpatiaLite"?

texts: When selecting fields to convert is there a way to select all within QGIS? I'm currently alternating space bar and a down arrow and this does not seem to be the most efficient method.


texts:
How can I create a circle (or "buffer") around a point anywhere on the globe in modern Open Layers with a radius specified in actual meters?
"Actual" meters meaning that a geodetic measurement from the center of the circle to one of its points would return about the same distance as was used to create the buffer. "Fake" meters like in EPSG:3857 are what I am trying to escape.
OL has a Circle constructor that takes a radius parameter which is used in units of the projection. As there is no global projection that allows us to create nice circular buffers with actual meters, this cannot be used.
Using Geodesic Measurements for Circle Radii? is for an hugely outdated version of Openlayers.


texts: Remove sections, defined by a line or polygon, from a raster in QGIS

texts:
I have created a raster layer of heights using the TIN interpolate function. I also have another raster layer with more accurate information on the height of particular areas defined by polygons/lines. I would like to remove these sections from the TIN generated raster and replace them with the information in the second raster.

I have tried using the clipper function, but this removes the wrong part (it retains the sections I want removed). Can anyone provide some guidance on how to achieve inserting values from a section of one raster into another?

QGIS 2.6.0 on Windows.


texts: "Back calculating" to quantify the "atmospheric effect" of water vapor or other gases in a "corrected" hyperspectral image using QUAC?

texts:
I know that when using FLAASH in ENVI, you can produce a water vapor map (caveats considered) from your atmospheric correction.

My question is if it is possible to "back-calculate" a water vapor map if I used the QUAC atmospheric correction? Possible for other gases (e.g. methane or CO2)?

I'm not so much interested in the accuracy of the output as much as I'm interested in the methodology. I would like to understand better, if only in theory, how FLAASH calculates its water map and if its possible for me to try my own calculation of atmospheric contributions for each pixel from the QUAC correction, which considers only 'in scene' radiance values.


texts: Problem in libtiff.dll when starting QGIS

texts: When I try to start QGIS 1.8.0 I get the error message "Can't find procedure start address TIFFMergerFieldInfo in libtiff.dll" and the program won't start. I have downloaded and reinstalled QGIS several times, but it does not help. What to do? Any suggestions? I have noticed that this problem is not unique to me, but I have not found a remedy.


texts: Split polygon into several smaller polygons in QGIS

texts:
I have a layer of polygons depicting response boundaries, I would like to break each of these polygons into a number of sub-polygons.  As this is for an exploratory project, the polygons themselves can be relatively random.

I have tried using Random Points and Voronoi Polygons to create a new mesh of smaller areas, however I need the sub-polygons to fit exactly within their parent polygon for each individual polygon (no sub-polygon should cross 2 parent polygons).

I have looked at Sub-dividing a polygon into smaller polygons, which seems relevant but I couldn't see how it applies, and Splitting large number of polygons at same time into smaller polygons? seems to do pretty much what I want, except my second layer is a polygon layer and the borders of the polygons do not match.

I am using QGIS, but also have access to SAGA GIS if that works better.


texts:
I want to compute the shortest paths travelling overland between the origin (Africa) and destination (North America).
My goal is to recreate human migration, and the shortest path should be across the Bering Strait.
First, I have converted a shapefile to a raster, imputing value 1 to land and a large value (999999999,000000) to water.
If my first step is correct, I do not have any no_data cell because I pre-initialize it with a value.

Second, I compute the cost using r.cost
function under GRASS.

Finally, I use r.drain to obtain the shortest path.

For some unknown reason to me, I must use the cost raster and the movement direction map ticking the The input raster map is a cost surface selector.

The result indicates that the shortest path consists in travelling via Iceland, thus crossing the Atlantic where each cell has a cost of 999999999 compared to the cost of 1 for travelling overland.

I think I cannot increase more the cost of water travelling.
How can I modify my workflow to force the Bearing route?


texts: Error code using 'Proximity (raster distance)' tool

texts:
I have been trying to use this function following 2 teaching methods and get the same problem.
I am trying to create a Boolean Raster layer from a previous Boolean Raster layer. The first boolean layer has 3 items from a collection of 5 that have been reclassified to group 3 items as 1 (true) with the remaining 2 as 0 (false).
I then used the proximity (raster distance ) tool following the instructions described by 2 different methods but get no output layer and the following error log:
GDAL command output:
C:\OSGEO4~1\bin\python3.exe: No module named gdal_proximity
Execution completed in 0.25 seconds
Results:
{'OUTPUT': 'C:/Users/user/AppData/Local/Temp/processing_0c64eca0c40547aba0fb1757949b2bbe/44b41e19a1ad45aeb316f25f65318b43/OUTPUT.tif'}
Loading resulting layers
The following layers were not correctly generated.C:/Users/user/AppData/Local/Temp/processing_0c64eca0c40547aba0fb1757949b2bbe/44b41e19a1ad45aeb316f25f65318b43/OUTPUT.tifYou can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.
Knowing me I am probably being stupid... but I have redone the exercise multiple times achieving the same result....
I have checked the bin and can't see the plugin/module with that name either.
FIXED..... So.... By uninstalling the LTR version I had (3.10.7) and following the procedure to install OSGeo4W so as to uninstall OSGeo4W, then physically removing OSGeo4W file location and shells, I then reinstalled a newer LTR version (3.10.8) and the problem appears to be fixed.


texts: Geographic coordinate systems on 2D screens

texts:
I know that a geographic coordinate systems are inherently three dimensional, while they are usually represented by two coordinates (lat and long) because the height is about the same across the surface.

These 3D coordinates are eventually projected to a 2D screen when we use a GIS software. My question is doesn't that mean they are practically projected? Isn't the only way a coordinate system is geographic, and not projected is that the map information is on a 3D sphere (like on a physical globe or Google Earth) and not on a 2D surface?


texts:
I've been trying to install add-ons wihtout success.  I keep getting an error message:

  ****IOError: [Errno socket error] [Errno 10060]**** Une tentative de
  connexion a échoué car le parti connecté n’a pas répondu
  convenablement au-delà d’une certaine durée ou une connexion
  établie a échoué car l’hôte de connexion n’a pas répondu
  (Mon Oct 28 22:16:44 2019) La commande s'est terminée (21 sec)

(Google Translation)

  ****IOError: [Errno socket error] [Errno 10060]**** A connection attempt failed because the connected party did not respond properly beyond a certain time, or an established connection failed because the login host did not respond (Mon Oct 28 22:16:44 2019) The order has ended (21 sec)

I've seen posts pertaining to this issue but haven't been able to fix it through editing the g.exetension.py file.  Is there any other way around this?


texts: Know of a source of historical 24-hour rainfall totals for US locations?

texts:
I'm looking to get historical information on the magnitude of rainfall events for a given location (Boston, Massachusetts USA). For example, to create a histogram of rainfall events by precipitation amount for a given location. So far, I've only found "last 3 days" or "last month" -type data sets from NOAA and others. I'd like many years worth of data if possible.

I'm not looking for a shapefile, so much as tabular-type data (database, Excel, CSV, etc).
Data doesn't need to be available for any arbitrary Lat/Lon, but could be for nearby POI, e.g. airports.


texts: How to change color of a layer in ArcGIS Javascript API by accessing the layer object?

texts:
I want to change layer color with every click of a button. So far I have been able to do lots working with the layer object alone, for example schools= this.map.getLayer("schools"). However, there are no methods for dealing with the symbology on the layer object.

I've seen documentation for changing symbology of graphics layers, however. But I'm not sure how to implement that from where I am now. How do I reference a graphics layer from my layer and access the fill methods? Or am I misunderstanding these objects to begin with?


texts:
I'm using dynamic text for owner names in my layout.  Problem is some are two words and some are a lot more.

When using a text box if the font is to big it just drops the extra.

Is there an easy way to simply have the font automatically resize to fit in the text box?  Something like the legend being able to resize/shrink to fit?


texts: What Algorithms to use for Levee breach Analysis?

texts:
I have a (relatively) High resolution DSM of the Country side, and want to figure out the extent of area that will be flooded if the Embankment/levee fails at a a given point.

Something like this: http://magic-gis.com/muscatineflood/index.html?appid=42cc730f8b20415d92ead71a89a22322&webmap=a90fb531890740c490ebeb25b2234876

I'm trying to do some research on this, and I'm drawing a blank, mostly because I don't know what this is called. Without proper keywords, A Google Search is extremely fruitless.

What is this kind of Analysis called? Are there any standard Algorithms/Models for this?

Are there any freely available tools for doing this kind of Analysis?


texts: Selecting records based on user input in ArcPy?

texts:
I am running Arc 10.3 in a versioned environment.  I am attempting to write a simple code that will select service records coming from Oracle CRM to a separate feature class containing only records that meet the SQL definition query that the user will enter. (Both feature classes, the Oracle CRM data and the new feature class CRM1 have the exact same field schema) The definition query will be a select records > than such and such date. I have the basics down but when I enter the SQL where clause it does not select any records.

Below is only the code for selecting out the records, the append part is not yet started as I can't get the select part to work first!

SQL Clause I am using :
"REPORTED_DATE" > date '2015-06-13 00:00:00'  (REPORTED_DATE is my field in both Oracle CRM and CRM1 that contains my date reported for the service requests)

Pic of my Model -

(source: gislab.net)

CODE BLOCK:

# Import system modules
import arcpy
from arcpy import env

# Set workspace
env.workspace = "Z:/Projects/sharedData/SMTEST.gdb/"

# Set local variables
in_features = "Z:/Projects/sharedData/SMTEST.gdb/CRM_ATTACH"
out_feature_class = "z:/Projects/sharedData/SMTEST.gdb/CRM_ATTACH1"
where_clause = raw_input()

# Execute Select
arcpy.Select_analysis(in_features, out_feature_class, where_clause)


texts: Edit features using drop down when using Edit widget in web appbuilder

texts: I have an app that displays the tree inventory of a university. The trees have categorical attributes like species, botanical name, etc and numerical attribues like height, diameter etc. When I use the Edit widget to edit the attributes of a tree, I have to type in the new species name. Instead, I would like to edit the categorical fields using a drop-down box that is populated from an excel file. Also, I would like to validate the numerical fields. Is there a way to do this?


texts:
ArcMap 10.2.2 for Desktop
License Type: Advanced
Python 2.7.5

Our SDE contains feature classes that use GlobalIDs for replication.  I am wanting to create a script that will run nightly backups on this SDE.  I have done plenty of scripts like this in the past but I haven't had to work around GlobalIDs before.

Testing:

If I simply go into catalog and highlight the feature classes that I need I can right click, copy, and paste into a new GDB and that will keep my GlobalIDs for all my features. But this is just a manual way of doing it that I know will work, but not what I need for automatic nightly runs.
If I do arcpy.FeatureClassToFeatureClass_conversion() this will keep the GlobalIDs field but generate new GlobalIDs for me.
If i do arcpy.CopyFeatures_management() this will remove the GlobalIDs field completely when moved into the new GDB.

I am hoping there is a way to mimic the manual 'right click copy/paste' method in python.


texts: Programmatical Setting of 'Store Relative Pathways' for multiple MXDs ArcGIS 10.1

texts:
I am looking to move large number of mxds to a different location on the same drive. Many mxds haven't got 'Store relative pathways' option checked therefore this operation will result in broken links to data. I've only started using ArcGIS 10.1 and it is now possible to programmatically update relative pathways for all mxds by using arcPy and setting relativePaths to be "True". I had a go at writing Python script which (in theory) loops through a top folder and searches for all mxds and then updates "Store relative pathways" setting to be on.
When I run this script in PyScripter, it fails to complete, I get the following error:
"exceptions.NameError: name 'mxd' is not defined

Can anyone help ?!?
Thanks in advance,
Magda

    import arcpy, os

folderPath = r"C:\My Folder Path\...\...\..."

for filename in os.listdir(folderPath):
    fullpath = os.path.join(folderPath, filename)
    if os.path.isfile(fullpath):
        basename, extension = os.path.splitext(fullpath)
        if extension.lower() == ".mxd":
            mxd = arcpy.mapping.MapDocument(fullpath)
            mxd.relativePaths = 'True'
        mxd.save()
del mxd


texts: Changing part of text element using ArcPy?

texts:
I work with multi MXD, and the title in all of them start with "The plan on the" and then i write the map's subject, for example:

"The plan on the land use", etc.
I need to change just the sentence from "The plan on the" into "The blueprint on the". All other words in the title will remain the same.
I try using this code:

import arcpy
from arcpy import env

env.workspace = r"D:\desktop\Project"
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"D:\desktop\Project\\" + mxdname)
        for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):
        if elm.text == 'The plan on the*':
            elm.text = 'The blueprint on the*'
            print 'elm.text'
    mxd.save()
del mxd

but it only print the MXD name:

>>>
airplane noise.mxd
airplane noise2 a3.mxd
antiquities.mxd
conflicts.mxd
project.mxd
>>>


texts: Counting unique species in polygon layer and buffer

texts:
Within QGIS I'm working on a data set that contains species at 100m resolution (eg OSGB 6 fig) and 10m or finer resolution (eg OSGB 8 or 10fig). The objective is to count the number of unique species within a polygon layer and associated buffer layer.

In acknowledging the resolution differences between the lower (6 fig) and higher (8-10fig) resolution, I have created a buffer around the polygons within which any 6fig points are also to be included (ie they fall outside the polygon layer but within the buffer).

While finding unique values for the buffer and the actual polygons is straightforward (count points in polygon), I'm not sure how to combine this for two different operations without inevitably double counting species that lie within both the polygon and the buffer layer. Similarly, counting all points within the buffer does not work as it also includes high resolution outside the polygon within adjacent polygons.

In essence, what I want is a way to combine the count points in polygon operation across two separate but linked features to find the total unique values within both the polygon and buffer layer without having duplicates of the same species included. One possible solution I have thought of is creating a presence/absence matrix for each species within the polygon/buffer layers then combining using field calculator, however I'm struggling to achieve this.


texts: Invalid QgsRasterLayer when attempting to load an image into QGIS - Mac (but not Windows)

texts:
From this post, Invalid QgsRasterLayer when attempting to load a Google Static Maps API image into QGIS this block of code can be used to load a Google Static Maps API image into QGIS (using the QGIS Python Console):

>>> reg = QgsMapLayerRegistry.instance()
>>> raster = QgsRasterLayer('http://maps.googleapis.com/maps/api/staticmap?center=51,0&zoom=10&size=2550x3300', 'test')
>>> raster.isValid()
>>> reg.addMapLayer(raster)

This works on Windows (tested versions QGIS 2.14.21 and QGIS 2.18.15). However, this does not appear to work on Mac (tested versions 2.18.9 and 2.14.20) or Linux (version 2.18.13). On Mac or Linux, and invalid raster object is returned.

This is not specific to the Google Static Maps API, an invalid raster object is returned with other http image urls.

Any ideas on how to get this working on Mac/Linux?


texts: QGIS geotrace languages from a given image of a map

texts:
I have a linguistic map of the Caucasus and I need to geotrace the languages that are on the map, preferably in QGIS.

Is there any way in which I can do it easily?

Here is the map:


texts: Increasing zoom level in 3D View QGIS Desktop?

texts:
I would like to increase the zoom level for my mosaic raster image in the 3D View of QGIS Desktop 3.4.4. I made a mosaic out of four raster images downloaded through the STRM-plugin for DEM data. Whenever I zoom in the 3D View, it seems as if the maximum size for me to zoom in is very small, so that I cannot increase it any further than this:

Does anyone know how to increase the zoom level for this image in the 3D View?

Could my problems stem from the image being a mosaic image?


texts: I have point in Google Earth Engine and I want to calculate the mean of the pixels around my point (point including). Now I do it manually (see photo), is there any function in GEE for this (any more professional solution)?


texts: A problem with Temporal, should it be classified as a bug?

texts:
I have previously used Temporal for animations with no problem but when I tried a new animation today with "Single field with date/time" the input below didn't show down arrow to select field, double L & R clicking made no difference..   A previously working animation has stopped working and shows same behaviour when temporal is selected.
The file in question had different columns for date and time.  So I modified file to first have date only column and second a merged date/time column.  Exactly the same result ie field selection arrow had no effect.  restarting QGIS between each attempt.  I suspect that QGIS temporal has become corrupted and will need to reload QGIS.
Unfortunately that is not the answer I removed/reinstalled QGIS but the original Temporal problem remains.  Any ideas?
The only way that this problem could arise is that
Temporal was modified on the lines of if field type not Date then ignore which would explain why a previously working animation failed.
I have tried the field calculator but the "Update field" is greyed out.

The only data required for an animation is lat, lon & date (or date/time) the other fields are not required.  So the solution is to use simple csv file in conjunction with csvt file.


texts:
I am trying to calculate canopy height. I have a RasterStack which has 8 rasters in it. The first 7 are DSM; the 8th is DTM. As you see in the function, I need to subtract every DSM with the same DTM. How can I specify x[[8]] this without saying 8?? or How can I use these rasters without stacking?? I want the code to be reproducible.
x in the function is rasterstack.
####Canopy Height Calculation####
CHM <- c()
chmCalc <- function(x) { #class(x) must be a "RasterStack"
  for (i in (1:(nlayers(x)-1))){
    res_chm <- x[[i]]-x[[8]] #DTM data needed to be set
    CHM <- c(CHM,res_chm)
  }
  names(CHM) <- paste0("CHM_H", 1:length(CHM))
  return(CHM)
} ```


texts: How to merge rasters with gdal_merge.py?

texts:
I'm trying to create a mosaic raster from 2 rasters. With rasterio works well but gdal_merge.py subprocess calling doesn`t work for me. MY ERROR IS :) : OSError: [WinError 193] %1 is not a valid Win32 application

import numpy as np
import matplotlib.pyplot as plt
import subprocess, glob
from osgeo import gdal

files_to_mosaic = glob.glob('C:/Users/DanielKovacs/Documents/GranCanaria/mosaic_L2A_T28RDS_A025431_20200505T115222_2020- 05-05_con/*_B04.tif')
files_to_mosaic

files_string = " ".join(files_to_mosaic)
print(files_string)

command = "gdal_merge.py -o mosaic2.tif -of gtiff " + files_string
output = subprocess.run(command)
output


texts: Single text label in dual font size in QGIS

texts:
I can create a point label with single size font (example)

However, I need to emphasize Japanese characters, make them larger in comparison to latin characters (example)

I tried some CASE WHEN expressions for Label -> Size option, but with no success.

Any hints (for QGIS version 2.18, alternatively v 3.x)?


texts: Using Spatial Join tool in ArcGIS Desktop?

texts:
I have a map with two layers: one of polygons and one of points (in coordinate pairs). To each point, I would like to match the attributes of the polygon in which it is located (or closest to). I am trying to use the Spatial Join tool. In the Target Features I have put in the points layer and in the Join Features I have put in the polygons layer.

The resulting layer has the coordinates, but unfortunately all attributes are <null>. Any idea why this is happening and how I can get the resulting layer to include the attributes as required?

Here are some screenshots to help understand the problem:
The layer of points and layer of polygons are:

And the resulting layer, after using the Spatial Join Tool as indicated above is:


texts: Can you access all tools in python/arcpy regardless of license level?

texts:
Is it possible to write a python script using arcpy that is for a tool not included in the current desktop license level?

For example, write something using the Identity Tool when you only have a standard license?


texts: arcpy Extract raster value to point: only extract values if point overlap the raster?

texts:
I have a set of 100 points and 100 rasters. Points are stored in one shapefile over large spatial extent and small, not-overlapping rasters are stored in gdb. Rasters are something like tiles. Together rasters cover the whole extent of points, but one by one, only one raster by one point can overlap. That is, many points are outside of particular rasters.

Here is an example: each point has its own one or two rasters, but rasters are not stored in mosaic, rather one by one:

I wonder how to Extract raster (multiple) values to points to for each raster to overlapping point?

Normally, I would firstly merge the individual rasters into new mosaic (Mosaic to new raster) and then simply Extract (Multi) Values to Points. NOw, I want to avoid the mosaic step due to memory issues (mosaic creates raster of 36 GB).

I was thining to something like this, using python and arcpy:

Pseudocode:

1. Read points shapefile
2. Read rasters into raster list
3. Loop through each row of shapefile (`arcpy.da.seach.coursor`)
4. Check if the shapefile row overlays one of rasters in raster lists?
    4.1. if yes - extract raster value to point, move to next row
    4.2. in no - move to next row

But I am not sure how to do this overlay check?? Or is there better approach?

Here is where I am now:

import arcpy, os
arcpy.CheckOutExtension("spatial")

# Define path to working directory
inWD = "C:/Users/extract"

# set working environment
arcpy.env.workspace = os.path.join(inWD, "analyzed.gdb")
arcpy.env.overwriteOutput = True

# Define input data
inPoints =  os.path.join(inWD,  "pts.shp")

# Read all raster in a list
rasters = arcpy.ListRasters("rst*", "GRID")

# Loop throught points row by row
with arcpy.da.SearchCursor(inPoints ,'OID@') as cursor:
    for row in cursor:
        # Create a query to subset rows one by one
        sql="""{0} = {1}""".format(arcpy.AddFieldDelimiters(inPoints ,
                                                            arcpy.Describe(inPoints ).OIDFieldName),row[0])
        print sql

        # Create new feature layer row by row
        outLyrName = "lyr_" + str(row[0])
        arcpy.MakeFeatureLayer_management(inPoints, outLyrName, sql)

        # Check if point overlays one of rasters:
        !! HOw to check this??

        # Process:
        arcpy.sa.ExtractMultiValuesToPoints(outLyrName, rasters, )

Simply extracting all points from all rasters in a list creates columns of empty data, so I think I can somehow do it one by one. I will highly appreciate advises how.


texts: Editing file geodatabase while other users read/view it gives ERROR 000224

texts:
Is there way to edit data within a file geodatabase while other user read/view data only. I deployed a file geodatabase and other users view the data but my process scripts can not edit the data and gives me ERROR 000224 cannot acquire a lock.
Is there a way to solve this problem?


texts: Missing 'Spec' Parameter using print.pdf on GeoServer Print Module

texts:
I have two configurations of Tomcat with the same version of GeoServer and the print module installed as below.

The old configuration is Tomcat 8.0.44/GeoServer 2.11.1
The new configuration is Tomcat 9.0.10/GeoServer 2.11.1

Using the old configuration I can POST to localhost:1080/geoserver/pdf/print.pdf in SoapUI with the spec file generated from the Print Module demo in the body and this returns a pdf.  An application has been written in C# based around this method of generating a pdf.

Using the new configuration (as Tomcat 8.0.x is now end of life, and also a combination of Tomcat 8.5.32/GeoServer 2.11.1) the same POST returns an error message:

<data contentType="text/plain" contentLength="55"><![CDATA[Error while
generating PDF:
Missing 'spec' parameter
]]></data>

If the POST is changed to a GET with ?spec={SPEC} then this works but for our actual application requirements a POST is the only way to ensure that all of the vector data gets sent to the print module.

Does anyone have experience of using Tomcat 8.5 or above with GeoServer's Print Module using a POST to print.pdf?

The default SPEC file that has been used in SOAPUI testing is

{
    "units":"degrees",
    "srs":"EPSG:4326",
    "layout":"A4 portrait",
    "dpi":75,
    "layers":[
    {
        "baseURL":"http://localhost:1080/geoserver/wms",
        "opacity":1,
        "singleTile":false,
        "type":"WMS",
        "layers":["topp:states"],
        "format":"image/jpeg",
        "styles":[""],
        "customParams":{}
    },
    {
        "type":"Vector",
        "styles":
        {
            "1":
            {
                "externalGraphic":"http://geoserver.org/img/geoserver-logo.png",
                "strokeColor":"red",
                "fillColor":"red",
                "fillOpacity":0.7,
                "strokeWidth":2,
                "pointRadius":12
            }
        },
        "styleProperty":"_gx_style",
        "geoJson":{
        "type":"FeatureCollection",
        "features":[{
            "type":"Feature",
            "id":"OpenLayers.Feature.Vector_52",
            "properties":{"_gx_style":1},
            "geometry":{
            "type":"Polygon",
            "coordinates":[[[-97,39],[-98,40],[-96,41],[-97,39]]]
        }
    },
    {
        "type":"Feature",
        "id":"OpenLayers.Feature.Vector_61",
        "properties":{"_gx_style":1},
        "geometry":{
            "type":"LineString",
            "coordinates":[[-97,40],[-98,39],[-99,38]]
        }
    },
    {
        "type":"Feature",
        "id":"OpenLayers.Feature.Vector_64",
        "properties":{"_gx_style":1},
        "geometry":{
            "type":"Point",
            "coordinates":[-98,38]
        }
    }
    ]},
    "name":"vector",
    "opacity":1
    }],
    "pages":[{
        "center":[-98.000000000002,40],
        "scale":4000000,
        "rotation":"0",
        "mapTitle":"A custom title",
        "comment":"A custom comment"
    }]
}


texts: Rotate bottom coordinates to be upside-down in print layout (QGIS)

texts:
Is it possible in QGIS (3.6) to somehow make all coordinates on the 4 sides of a map in a print layout point towards the map?

I found the option for the 2 vertical sides. If the map is rotated by 90° because of space and annotation issues, it would be nicer to have the bottom coordinates upside-down, pointing towards the map.

Then it wouldn't matter so much how the print out is rotated.

I thought of adding two grids, one for the bottom labels only. I can't find any option to rotate the labels independently of the orientation setting though.


texts:
I'm trying to publish my VB.net Application, it's working fine on my machine, but it's showing this error on my friend machine:

  unable to load dll 'arcgisversion.dll': The specified module could not
  be found.

I've checked platform and target cpu, they both are 86x.

What I have tried already:
 - releasing the project, which adds alot of dlls, the error remains
 - copying arcgisversion .dll from my pc to the target pc, and applying this command as Administrator :regsvr32 /i ArcGISVersion.dll, an error message shows up: "the module arcgisversion was loaded but the entry-point dllinstall was not found"


texts:
I have a polyline layer that has been joined to a table and I need to split it at certain distances along the line.

The attribute table contains the distances at which I would like to split the lines. The attribute table looks like this:

I need to split the line a two spots, ChStaL and ChStaR, which are both distances from the starting node.

I'd really like a python type solution if possible as I have over 30 of these layers.

I've been cruising around ESRI help docs and ArcToolbox but, haven't found anything that will allow me to do this en-mass. Almost everything out there requires me to do manually editing.

Any help would be most appreciated.


texts: Filter MODIS data for particular sensor zenith angle using earth engine

texts:
I need surface reflectance data with a sensor zenith angle <45 degrees. My current code for
green ndvi is: I have selected the related surface reflectances and the sensor zenith angle layer. But I don't know how to filter for observations with the view angle<45. Can someone help with this please.

 var A = ee.ImageCollection("MODIS/006/MYD09GA"),

 var NIR = A.select('sur_refl_b02')
 var Green = A.select ('sur_refl_b04')
 var sza = A.select('SensorZenith')

var addGNDVI = function(image) {
var gndvi = ee.Image(image).normalizedDifference(['sur_refl_b02', 'sur_refl_b04']).rename('G');
return ee.Image(image).addBands(gndvi);
}

//Map the GREENNDVI over collection A
var AGNDVI = A.map(addGNDVI)
var AGNDVIdated = AGNDVI.filterDate('2012-06-01','2014-09-30')


texts: Displaying GeometryCollection with one geometry in QGIS fails

texts:
I have a problem with displaying simple geometrycollection with one multipolygon as a vector layer in qgis. When I import the geometry in qgis I cannot zoom to layer nor see labels and style.
Here is the link to data I want to show GeometryCollection

I managed however to show the data in geojson.io and tried to convert to different formats (shapefile,csv) but converted files were empty.

What would be the issue?


texts: Adding Z value from a Point field to a LineString using pyqgis

texts:
I Have two layer in Shapefile format:

a line layer made from a set of point which is a polyline 3D with z value to 0.
a point layer which contains elevation in a field (field_4)

Each point is located on a vertex of the line. See below:

What I want to do is to had the Z value containing in the field_4 of the point layer to the coresponding vertex.

I found the setZAt method from QgsLineStringV2 which seems to exactly fit my need. So I wrote this code snippet :

line_layer  = QgsMapLayerRegistry.instance().mapLayersByName('line_truck')[0]
point = QgsMapLayerRegistry.instance().mapLayersByName('RT03')[0]
line_layer.startEditing()

for feat in line_layer.getFeatures():
    for i, elem in enumerate(point.getFeatures()):
        #get the point lying on the vertex
        if str(feat.geometry().vertexAt(i).x())+str(feat.geometry().vertexAt(i).y()) == str(elem.geometry().asPoint().x())+str(elem.geometry().asPoint().y()):
            print i , elem['field_4']
            #update z value but nothing happen...
            #feat.geometry().geometry is a QgsLineStringV2
            feat.geometry().geometry().setZAt(i, elem['field_4'])

When I run it, the print statement show me the right values but the z value isn't updated in my feature...

Do I miss some steps or make a mistake in my logic??


texts: Creating spatially projected polygon grid with ArcMap?

texts:
i need to create a grid with 1x1 square kilometres cells in ArcMap. The operation located in  Data Management Tools -> Feature Class -> Create Net is almost what i need but in the parameters, i can't define a projection the raster grid ist based on. The raster grid must be in the Projection UTM 32N, not in an "undefined" projection.

Is there any operation/function in ArcMap i can use to create a spatially projected polygon grid?


texts: Reprojecting unequal number of latitude and longitude: pyproj.exceptions.ProjError: x, y, z, and time must be same size

texts:
I have a rectangular matrix of geographic climate data with a 1d row of latitude coordinates (size=834) and a 1d column with longitude coordinates (size=1115). They are both in unit of m (EPSG:32630) and I would like to transform to lat/lon. When I try to do:
    from pyproj import Proj, transform

    inProj = Proj('epsg:32630')
    outProj = Proj('epsg:4326')
    x1 = precip.lon.values
    x2 = precip.lat.values
    x2, y2 = transform(inProj, outProj, precip.lon.values, precip.lat.values)
    print (x2, y2)

I get the error:
Traceback (most recent call last):
  File "<input>", line 4, in <module>
  File "<conda_env>/lib/python3.7/site-packages/pyproj/transformer.py", line 759, in transform
    ).transform(xx=x, yy=y, zz=z, tt=tt, radians=radians, errcheck=errcheck)
  File "<conda_env>/lib/python3.7/site-packages/pyproj/transformer.py", line 446, in transform
    errcheck=errcheck,
  File "pyproj/_transformer.pyx", line 459, in pyproj._transformer._Transformer._transform
pyproj.exceptions.ProjError: x, y, z, and time must be same size

Does anyone know how I can transform from PCS to GCS with Python without having equal amount of latitude and longitude values?


texts:
I'm a bit of newbie at GIS so be gentle with me.

In QGIS I've set up a WMS from the British Geological Survey:

http://maps.bgs.ac.uk/ArcGIS/services/BGS_Detailed_Geology/MapServer/WMSServer?

Now the layers are rendering at 1:25000 but at 1:50000 and above the screen goes blank.  According to the BGS website they should be visible between 1:75000 and 1:10000. Does anyone have any suggestions as to what I'm doing wrong?


texts:
I installed Orfeo using OSGeo4W and have set paths in QGIS (Wein) at

Processing -> Options & Configuration -> Providers -> Orfeo Toolbox

to:

Applications: C:\OSGeo4W\apps\orfeotoolbox\applications

Command line: C:\OSGeo4W\bin

And the contents of these folders include the required files.

When I try to run a command I get an immediate error:

  otbApplicationLauncherCommandLine.exe has stopped working

The log shows no errors and appears to have executed.

I wonder if I'm missing an installation step. I can't see Orfeo in the plugins menu for example. I'm running Windows 8.1, and used 64 bit OSGeo4W installer.


texts: Image classification from "Imagery" basemap ArcGIS 10.2

texts:
I'm trying to reclassify the coast of France and Spain into rocky or sandy categories. So far, I've been drawing polygons around the different patches, but want to try using image classification to speed up the process.  Is there a way I can use the World_Imagery basemap as a normal raster which I can then use in image classification?

There was a similar question where the only answer was to "Export Map", but this would take as long as doing all the polygons by hand.


texts:
In my code I'm putting on an ol.interaction.Select interaction. I need to set the style and this is easy.

   new ol.interaction.Select({
            filter : function(feature, layer) {
                return me.inFilterSelect(feature, layer);
            },
            wrapX : false,
            style : stylePoligon
        });

I now have the need to set two styles, but one must be a function style (the styles depends on some attributes attached on the feature), and the other could be a function too or a simple new ol.style.style({...}).

How can I set such a thing, multiple styles of different types (one function and one ol.style.style) or multiple functions style? I'm using openlayer 4.0.1.


texts: What are some high-profile websites that use OpenLayers and PostGIS?

texts:
I think there are many convincing arguments for using PostGIS but I am interested in showing people examples of websites that use PostGIS and OpenLayers.

Some examples that I have seen are:

foursquare
CartoDB
Vizzuality (many examples)
WalkScore (as far as I can tell)

Are there other examples that show good examples of this combination? I'm not sure if any of the sites that I listed are using OpenLayers as a central part of their functionality but I believe that they all use PostGIS.


texts: Summary Statistics with multiple case fields not working

texts:
While trying to calculate summary statistics in a table in a GDB, i get the following error:

ERROR 999999: Error executing function.
The table was not found. [OutputStats]
Failed to execute (Summary Statistics).

The "OutputStats" table is the output of the Summary Statistics tool.

The input table table is in the following format:

Address       | Period | Consumed
1 Example str | 1/2016 |    18
1 Example str | 2/2016 |    17
2 Example str | 1/2016 |    19
2 Example str | 2/2016 |    23

The table has 300k lines but i have worked with larger ones without problem.

I use "Consumed" as the statistics field with "sum" as Type with "Period" and "Address" as the case fields.

Note that "address" and "period" are text fields and "consumed" is Numeric.
When i use only one of the above case fields the tool works ok.
I have tried both ArcGIS 10.5 and ArcGIS Pro 1.4 with no success.
I am calling the tool in Model Builder but i also tried it from the toolbox as a standalone tool.
In Model Builder i also tried a Table View as input without success.

I think the problem is with the input table but i can not see why.


texts: Export Joined table an Feature layer based on a query

texts:
I have a table and a layer joined fine, I am trying to export the end result into a file geodatabase, the export is fine if I specify specific name, other than that it exports the original layer data.
Here is the code I am using:

public void JoinLT(string DataPath,string Sheet, int lyrIndex, string JoinName, string TableFieldName,string FCFieldName, string FCName, string TableName)

    //Export the lakes from the GNE_Layer, to an already defined file geodatabase
    {
        try
        {
            // Part 1: Define the Environoment
            IMxDocument mxDoc = (IMxDocument)ArcMap.Application.Document;
            IMap map = mxDoc.FocusMap;
            ILayer layer = map.Layer[lyrIndex]; // GNE_Clip
            IFeatureLayer featureLayer = (IFeatureLayer)layer;
            IFeatureClass featureClass = featureLayer.FeatureClass;

            IWorkspaceFactory WsF = new ExcelWorkspaceFactory();
            //Open the Excel sheet
            IFeatureWorkspace fWs = (IFeatureWorkspace)WsF.OpenFromFile(DataPath, 0);

            //Open table "tableName"
            ITable table = fWs.OpenTable(Sheet); //Dynamic
            //ITable table = fWs.OpenTable("TroutLakes$");

            //Add StandaloneTable
            IStandaloneTableCollection ISTC = (IStandaloneTableCollection)map;
            IStandaloneTable SAT = new StandaloneTable();

            SAT.Table = table;
            ISTC.AddStandaloneTable(SAT);

            //Refresh Activeview
            mxDoc.ActiveView.Refresh();
            mxDoc.UpdateContents();

            // Part 3: Join the table to the feature class
            // Create a memory relationship class
            IMemoryRelationshipClassFactory memRelFact;
            IRelationshipClass relClass;
            memRelFact = new MemoryRelationshipClassFactory();

            relClass = memRelFact.Open(JoinName, featureClass,
               TableFieldName, (IObjectClass)table, FCFieldName, FCName, TableName,
               esriRelCardinality.esriRelCardinalityOneToMany);

            // Perform a join
             IDisplayRelationshipClass dispRC = (IDisplayRelationshipClass)featureLayer;
            dispRC.DisplayRelationshipClass(relClass, esriJoinType.esriLeftInnerJoin);

            IDisplayTable dispTable = (IDisplayTable)dispRC;
            IDatasetName srcName = (dispTable as IDataset).FullName as IDatasetName;
            IDatasetName destName = new FeatureClassNameClass() as IDatasetName;

            string LTOutput = Interaction.InputBox("Please enter Feature class name");
            destName.Name = LTOutput;

            IWorkspaceName destWsName = new WorkspaceNameClass();
            destWsName.PathName = "Data\\trial.gdb";
            destWsName.WorkspaceFactoryProgID = "esriDataSourcesFile.FileGDBWorkspaceFactory";

            destName.WorkspaceName = destWsName;

            IFeatureClassName destFcName = destName as IFeatureClassName;
            destFcName.FeatureType = esriFeatureType.esriFTSimple;
            destFcName.ShapeFieldName = "Shape";
            destFcName.ShapeType = featureLayer.FeatureClass.ShapeType;

            var queryFilter = new QueryFilterClass();
            queryFilter.SubFields = "*";
            queryFilter.WhereClause = "(GNE_Clip.Lake_Town = \"TroutLakes.Lake_Town\")";
            ISelectionSet selnSet = null;
            IGeometryDef geomDef = null;

            IExportOperation exOp = new ExportOperation();
            exOp.ExportFeatureClass((IDatasetName)srcName, queryFilter, selnSet, geomDef, destFcName, 0);

            mxDoc.ActiveView.Refresh();
            mxDoc.UpdateContents();

        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.Message, ex.Source);
        }
    }


texts:
I have a large osm dataset. When I execute a pgr_ksp query, it takes +30 minutes

Example:

SELECT seq , id1 , id2 , cost
 FROM pgr_ksp(
   'SELECT id, source, target, cost FROM myways',
   524356, 1210854, 2, false);

My virtual server has 4 processors and 3GB RAM. My ways table has 2290470 rows.

The other functions work normally but pgr_ksp does not.

Does anyone know what is the problem?


texts: Can I move selected features in ArcMap with arrow keys?

texts: I would like to be able to move features by small (preferably predefined) increments using the keyboard arrow keys?  Can this be done?


texts: Updating field in attribute table of multiple shapefiles in ArcGIS 10.1

texts:
I have almost 100 shapefiles (all having polyline features only) in a single folder with each filename "0","1"..."99" . I need to recalculate the length and another field "slope" which I previously calculated in field calculator:
slope =
math.degrees(math.atan(100/ !Length! ))".

Now I want to recalculate both fields in each shapefile.
PS: These 100 shapefiles were previously one shapefile whose length and slope has been calculated. But there after I clipped it using features of another shapefile resulting into 100 shapefiles. But the length and slope value was not get upated and it is showing the old values of the single shape file. Now I want to calculate slope and length again for each shapefile which is time consuming so I hope something can be done in a Model Builder or a Python script which could ease my computation.


texts: Cannot get areas for polygon layers in ArcMap?

texts:
For an ArcGIS project opened in ArcMap 10, we have added polygon layers representing different geographical features.

After all the polygons were created, we opened the attribute tables for those layers, and to our shock/dismay, the cells in the "shape_area" column all showed zero! (The "shape_length" field was correctly filled.) We also tried the Measure tool in the toolbar, but it could only measure distances, and the Area button was greyed out.

We checked to make sure that all layers were in the correct geographical coordinate system (WGS 1984).

What did we do wrong? Is there any way to obtain the areas of the polygons we created?

As Brad suggested, I tried to right click on the "SHAPE_Area" column (both on the column label and in the column), but the "Calculate geometry" (and "Field calculator") options were also greyed out.


texts:
I want to use gdal_retile. I have gdal installed after installing it from OSGeo4W by my bin folder is missing gdal_retile. How can I install it? I am working on a Windows machine. Is it expected not to be part of GDAL from OSGeo4W?
Edit: I can see my question has been closed due to duplication.
What I am confused about is how I am suppose to install and use gdal_retile. I found an OLD installation from OSGeo4W which include bat files so I can run it as a command line tool. But that does not seem to be the case anymore? gdal_retile itself is a python script? But then I must run it using python not a command line tool? I am completely stumped in what direction I should go to get it running. If I want to use gdal_retile I assume it is not enough to install just gdal?
You can install GDAL in numerous ways
https://gdal.org/download.html
Using the OSGeo4W installer
Or binaries from https://www.gisinternals.com/release.php
I am not a computer scientist and there are too many options for me to understand the differences between them and the requirements to get going. I can post a new question if that make more sense.


texts: how i calculate land surface temperature using landsat 8 in Envi 5.1

texts: I am doing research about land surface temperature using Landsat 8 image. However, I have a problem with them. I have read document from Landsat.org, but i can not do it in Envi. Could you help me step-by-step to calculate land surface temperature using landsat 8 TIRS in Envi 5.1.


texts: Arcpy UpdateCursor 'SHAPE@X', 'SHAPE@Y' in decimal degrees

texts:
I have the code below that works, however I want the units to be in decimal degrees. How do I force this?

fields2 = ['cent_lat', 'cent_long', 'SHAPE@X', 'SHAPE@Y']

with arcpy.da.UpdateCursor(fc,fields2) as cursor:
    for row in cursor:
        row[0] = row[2]
        row[1] = row[3]

        cursor.updateRow(row)


texts: I have been creating some QGIS plugins lately and configuring the interface of QGIS the way I like it in a programme. I would now like to hand over the system to other people: Is it possible to create an exe-file or something similar which automatically installs QGIS with the correct settings and plugins on another computer?


texts: Geodesic distance from raster in python?

texts:
I am carrying out some analysis in python (using Anaconda), I have a binary raster and want to perform a geodesic distance calculation from all areas with a value of 1- I can do this within QGIS using the proximity raster tool. In python I have been using gdal, however the gdal_proximity function that QGIS uses does not seem to be available, I have been using the gdal documentation https://gdal.org/python/ and cannot see it there either.

Is this just unavailable when using gdal in python and, if so, is there an alternative?


texts: Stratified Random Sampling with large dataset in ArcGIS?

texts:
I want to create a Stratified Random Sample with a land cover classification map. The purpose for the Stratrified Random Sample is to produce an Accuracy Assessment. I am using ArcGIS 10.2 with access to Spatial Analyst and all other tools.

I have six land cover types and I want 300 random points for three of the classes I am interested in.My first step was converting my raster to polygon. The result gives me over 3 million polygons to work with.

I have tried the Dissolve tool and making multi-part polygons, but some of my classes are too large, which results in an error message "Topology Invalid (Out of Memory).

I have tried using the Merge in the Editor toolbar. Again, due to my large number of polygons in some of my classes, it freezes and crashes ArcMap.

I would like to do the stratified random sample on a raster, rather than converting it to polygon form. Is this possible?

Are there any other tools/programs that I should use that may help combat my large dataset issue? I have consistently run into problems because I simply have too many polygons per land cover type.


texts: Specifying layer and BBOX in URL when opening GeoExplorer?

texts:
I would like to specify the state that the GeoExplorer map window and which layers are turned on when opening it up by a URL.

I have started to go through the code but I feel this must be quite a common user requirement.

How can I do this?


texts: Coloring thousands of polygons so each is distinct from its neighbors?

texts:
I am using ArcGIS 10, and I have layers of the 2010 Census blocks, block groups, and tracts.  There are thousands of polygons in these layers.

I want to color them so that each is distinct from its neighbors, but all I can figure out is how to give each one a distinct color.  The color of each polygon does not really matter; I just want to easily distinguish each from its neighbor by color, rather than with a border.

Or is that something you only see in an elementary school map of the 50 United States?


texts:
I have downloaded a full set of text DTM's for Northern Ireland from OSNI Open Data, but when I load them into QGIS, some of them do not appear. When I try to load the missing files, I get an error message:

"Ungridded dataset: At line 38676, too many stepY values" (this is from one of the files, so the line number may differ for each offending file).

The URL to download one of the sets is:
https://opendatani.blob.core.windows.net/lpsopendata/10M-DTM-Sheets1-50.zip

In this dataset, you will find that 'sheets' 2, 6 and 16 are the problem files.

Is anyone able to shed any light on what may be causing the problem?


texts: Unable to open datasource with ogr2ogr.main()

texts:
I try to make ogr2ogr work in Python the way suggested by Elevine's answer elsewhere, but I'm facing bugs. Most notable error is that while I try to reproject shapefiles, I get the following error:

Unable to open datasource `lesy_w.shp' with the following drivers.

A list of formats follow, including "ESRI SHapefile". My command looks like this:

ogr2ogr.main(["", "-f","ESRI Shapefile", "-s_srs","epsg:4326", "-t_srs","epsg:32633", "lesy_u.shp","lesy_w.shp"])

The corresponding ogr2ogr command (ogr2ogr -f ESRI Shapefile -t_srs epsg:32633 lesy_u.shp lesy_w.shp) works well in the command line. It's a simple line geometry, no "geometry collection". The input shapefile (lesy_w.shp) is loaded as a QGIS layer normally. I tested this with other layers and the problem is still there.

I get the same error also while trying other operations, like converting to KML without reprojection. The same when I try to convert from kml to shapefile.

Is filing this as a bug in ogr2ogr.main() my only option, or is there any solution?


texts: QgsVectorFileWriter returns false output

texts:
I am trying to make a script that reads the content of a .dbf-file and exports it as a shapefile. Equal to the manual pattern "Save as.." when rightclicking the imported layer in the Qgis Desktop.

I've read the Cookbook for developing in pyqgis and tried to troubleshoot with the documentation of QgsVectorFileWriter() but it seems to be written for C++ implementation and not for Python.

 _VctLay = QgsVectorLayer("path-to-file\file_imp.dbf", "Import_Area", "delimitedtext")
 QgsVectorFileWriter.writeAsVectorFormat(_VctLay ,"path-to-file\file_exp","CP1252",None,"ESRI Shapefile")

I'd like it to be a Shapefile but it is only the .dbf-file i have to use. Moreover I want the CRS to be "EPSG:4095" but it makes an error if i type it in like that. I've tried to have the "bool onlySelected = True" and False" - but none is adding content to the exported table.

I am getting an output. The whole package of files that follows with a shapefile. But they are about all empty. It only creates the files themselves.

Is there some flag, some process i am missing out on this? I wonder if i could have the first step. _VctLay = QgsVectorLayer("path-to-file\file_imp.dbf", "Import_Area", "delimitedtext") to be followed by an print-commands that clarifies the content of the imported .dbf. But exactly how, i don't know yet.


texts: Spatial query in PostgreSQL to find a route to a specific point?

texts:
I drew 2 shapefiles - point.shp and line.shp. Line.shp has some lines intersecting points in point.shp. I converted the shapefile to PostgreSQL, but how can I make a spatial query to select a line from one point to another point ?

Thanks all.


texts: Facing a problem when re-projecting a file?

texts:
I am trying to re-project from Equi-rectangular. The file is provided as HDF5.
I wonder why half of the map is not projected correctly.

gdal_translate -of "ENVI" -gcp 0 0 0 90 -gcp 1440 0 360 90 -gcp 0 720 0 -90 -
gcp 1440 7 20 360 -90 -a_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

here how the results look like:
![enter image description here][1]

For info:


texts: How to set the active tool in the composerview?

texts:
This script works fine for me, but when setting the active tool, its active but not marked as active in the gui (still the default tool is highlighted):

from qgis.PyQt.QtXml import QDomDocument
from qgis import gui

scale = 10000
templateFile='c:/1temp/mylayout.qpt'
activetool = 15# 'MoveItemContent'

composerview = iface.createNewComposer()
#load template
f= file(templateFile, 'r')
templateContent = f.read()
f.close()
doc=QDomDocument()
doc.setContent(templateContent)
composerview.composition().loadFromTemplate(doc)
#show composer
composerview.show()
#setup map
map = composerview.composition().getComposerMapById(0)#QgsComposerMap
#extent
map.zoomToExtent(iface.mapCanvas().extent())
#scale
if scale >0:
    map.setNewScale(int(scale))
#set activetool
if activetool is not None and activetool<>'':
    composerview.setCurrentTool(gui.QgsComposerView.Tool(activetool))
#finally
composerview.composition().refreshItems()


texts: Vectorizing a raster zone map without the lines turning into squares QGIS

texts: Raster zone lines turn into squares after vectorizing my raster because of pixels between color zones QGIS


texts: Add AutoCAD .dwg annotation file to Arcmap using python (Same as +ADD Data button)

texts:
I would like to add an AutoCad .dwg file to Arcmap using python the same way that the +ADD Data button does.

I can add an Annotation layer with python but only as a feature layer not as an annotation layer... even thought in the properties it says it's a CAD Annotation Feature Class

targetGroupLayer = arcpy.mapping.ListLayers(mxd, "ACAD", df)[0]
addLayer = arcpy.mapping.Layer(PlanPath)
arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM")

The files share the same properties, but are different types (Feature, Annotation)...

NOTE: I realize in my screencaps the files have different R1, R2 numbers ignore this.

The documentation for layer says:

  "There are essentially three categories of layers in a map document: feature layers, group layers, and raster layers. There are a few specialized layers and datasets that don't fall into one of these three categories: annotation subclasses Group layers and other sublayers (e.g., annotation classes) are treated just like ordinary layers."

Is there any way to add it as an Annotation layer? The same way the +Add Data button works?

EDIT:

This is what I tried doing.

targetGroupLayer = arcpy.mapping.ListLayers(mxd, "ACAD", df)[0]
arcpy.MakeFeatureLayer_management(PlanPath,"CAD")

for lyr in arcpy.mapping.ListLayers(mxd,"CAD",df):
    lyr2 = lyr
    #Include the Add Layer to Group code as part of the For Loop
    arcpy.mapping.AddLayerToGroup(df,targetGroupLayer,lyr2)

From what I could tell from your code (Get Spatial). The annotation layer is added to a temp feature layer. You then search for that temp feature layer "cadname" and then add it to the target grouplayer.I tried to recreate this but using and for loop to loop through each layer and find the one called "CAD"

@GetSpatial I ran the new code and looks like it made the feature layer yet did not correctly add it to the map.

In the report below is the "CAD" after "\Annotation" part of the file name it is trying to find, or is it simply the name of the Layer and has no effect on the file path?


texts:
I want to write scripts for Erdas in Python.

I downloaded some Python examples here but I could not make much.

Where can I get the full Help material for Erdas python scripting?


texts: Can anyone give me a hand to install the Mapbox studio in Ubuntu 12.04? I have downloaded the zip file for Linux but cannot find any installation instructions there.


texts: Function that applies another function to nearest neighbors of a specified location?

texts:
I'm sure this is quite basic, but I just couldn't find anything appropriate: I am looking for a function that applies another function to variables of a SpatialPointsDataFrame, but only for locations that are the nearest neighbors of a specified location.

Alternatively, how can I subset a SpatialPointsDataFrame by distance to a given location?

For example, in the meuse data set (gstat package) I would like to calculate mean, variance and quantiles of zinc concentration for the 20 closest locations that are at most 0.1km away from location x=179997, y=331662.

I found the function applynbd() in the spatstat package, but it applies the given function to all points in the data frame insteat of a single one. Moreover, only a number of nearest neighbors OR a distance (max. radius) can be specified instead of both.


texts: I am developing a set of plug-ins for QGis (2.14). Every time I need to try it out I must quit and relaunch QGis. Is there any way to force QGis to reload the plug-ins without restarting it?


texts: Labeling all features of shapefile with PyQGIS?

texts:
I have a problem when trying to label all features of a shapefile with a Python routine.

If I open my shapefile in QGIS "by hand", and select a .qml with pre-defined label selection, the labels are correctly displayed. If I load the .qml file via a Python routine (see line below), all style features (colors etc.) are there, EXCEPT for the labels. No idea why.

myLayer.loadNamedStyle("someStyle.qml")

Since the line above did not give me the desired labels, I tried to create a more extensive routine, see below.

myProvider = myLayer.dataProvider()
myLayer.startEditing()
for myFeat in myLayer.getFeatures():
   #Some calculations in here...
   #The variable "Some_Value" is what I'd like to use as label

   Index_myLabelField=myLayer.fieldNameIndex("myLabelField")
   Value_myField={myFeat.id():{Index_myLabelField:QVariant(Some_Value)}}
   myProvider.changeAttributeValues(Value_myField)
   myLayer.commitChanges()

   myLayer.label().setLabelField(0,Index_myLabelField)
   myLayer.label().labelField(0)
   print myLayer.label().fieldValue(0,myFeat)   #Returns the desired label string

   #Create map renderer
   mapRenderer = QgsMapRenderer()
   rect = myLayer.extent()
   mapRenderer.setExtent(rect)
   mapRenderer.setLayerSet(QStringList(myLayer.id()))
   myRenderContext = mapRenderer.rendererContext()

   myLayer.select(myFeat.id())
   print myLayer.isValid()  #Returns true
   myLayer.label().renderLabel(myRenderContext,myFeat,1)   #Error: Segmentation fault (core dumped)

Until the last line, where I try to render the label, everything seems to work fine.

What is going wrong in the last line?

There are two things I am unsure about:

In the QGIS API, the function to tell QGIS the field I want to use as label is

setLabelField (int attr, int fieldIndex)

What is the "attr" in this context?

What should I put here?

I used "0" so far...

In the renderLabel() function, there is a "bool selected" needed in third position of the call. So what I am doing is selecting the feature (myFeat), and telling the function a "1" here.

Is that correct?

Rendering and Labeling Shapefile with PyQGIS does not fully help me, as the OP there does not seem to label all features of his shapefile. I am able to add one label to the shapefile (kind of like a title) with the described method there, but I am looking to have a label next to every feature (e.g. like street names on a map).

If you know why the qml-loading in the beginning of my post does not work can you let me know?


texts:
The following script creates a new version in SDE but generates the following error below. Am I missing an environment?

# Name: SDE_CreateVersion
# Description: Creates a new SDE version

# Import system modules
import os
import arcpy

# Set local variables
sdeConnection = "Database Connections\Jordan.sde"
parentVersion = "sde.DEFAULT"
versionName = os.environ.get( "USERNAME" )

for version in arcpy.da.ListVersions(sdeConnection):
    if version.name.split(".")[0] == versionName:
        break
    else:
        # Execute CreateVersion
        arcpy.CreateVersion_management(sdeConnection, parentVersion, versionName, "PUBLIC")

This is the error produced:

Traceback (most recent call last):
  File "C:\Users\jnmiller\Desktop\Create Version.py", line 20, in <module>
    arcpy.CreateVersion_management(sdeConnection, parentVersion, versionName, "PUBLIC")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\management.py", line 18856, in CreateVersion
    raise e
ExecuteError: ERROR 001148: Cannot create a version with the specified name.
Failed to execute (CreateVersion).