texts: Extract coordinates from polyline in QGIS

texts:
I have a shapefile with 64 lines with the metadata:

General:
Storage type of this layer: ESRI Shapefile
Source for this layer: C:/Users/toke.nielsen/Desktop/lines.shp
Geometry type of the features in this layer: Line
The number of features in this layer: 64
Editing capabilities of this layer: Add Features, Delete Features, Change Attribute Values, Add Attributes, Delete Attributes, Create Spatial Index, Fast Access to Features at ID, Change Geometries
Extents:
In layer spatial reference system units : xMin,yMin 563677.42,6223900.13 : xMax,yMax 568136.66,6227899.44
Layer Spatial Reference System:
+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

How do I export the coordinates from these lines to a .csv file? I have tried to write GEOMETRY=AS_XY in the layer options when creating the .csv and I have tried the vector->geometry tools->export/add geometry columns, but then I only get the length of each line.


texts: Efficient method for determining directional change in landcover polygons between two time points?

texts: I have polygons of landcover classes in ArcMap (dense forest, sparse forest, wetland, water, rock) for a landscape from two different time points. I want to determine the % or total area that changes from one landcover class to the others between the two time periods. I know that I can likely do this all feature-by-feature using the "Intersect" tool, but is there a tool or extension that would streamline this process?


texts: ModuleNotFoundError: No module named 'osgeo' in virtualenvironment

texts:
As I am working in ubuntu. This works if I deactivate the virtual environment but after activating the environment I got the error saying no module named osgeo. Any help would be appreciated to install it.

from osgeo import ogr, osr


texts:
I have made an ArcPy add-in that assists a user by adding layers from a nominated fGDB to their map based on a nominated extent. I wish to add a button that changes the data frame properties to clip all the TOC layers to shape (the extent polygon). I can't find an ArcGIS Desktop function to do this. There is however a production mapping (PLTS Extension) function that might suit, but we don't wish to purchase a "solution" for a 9 click standard desktop process.

ArPyProduction Clip To Data-frame:
http://resources.arcgis.com/en/help/main/10.1/index.html#//010300000306000000

Any ideas?

class ButtonClass3(object):
    """Implementation for Version1_Addin_addin.button_1 (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
    print "Limiting display of layers to your extent"

    mxd = arcpy.mapping.MapDocument('current')
    df = mxd.activeDataFrame
    global ExtentFC
    ExtentFC = "C:\\Data\\ArcGIS\\SelectData.gdb\\ExtentFC"
    global ExtentFCLayer
    ExtentFCLayer ="Current_Extent"
    global ExtentLyrPath
    ExtentLyrPath = "C:\\Data\\ArcGIS\\ExtentLayerFile.lyr"

    ### Carto
    arcpy.MakeFeatureLayer_management(ExtentFC, ExtentFCLayer"", "", "")
    FeatureClassLayerMapping = arcpy.mapping.Layer(ExtentFCLayer)
    ExtentLyr = arcpy.mapping.Layer(ExtentLyrPath)
    arcpy.mapping.UpdateLayer(df, FeatureClassLayerMapping, ExtentLyr, True)

    ### Select Extent & Zoom
    arcpy.SelectLayerByAttribute_management(ExtentFCLayer,"NEW_SELECTION","\"OBJECTID\" >0 ")
    df.zoomToSelectedFeatures()
    arcpy.RefreshActiveView()

    ### arcpyproduction EXAMPLE A:
    arcpyproduction.mapping.ClipDataFrameToGeometry (data_frame, clip_geometry, {exclusion_layers})

    ### arcpyproduction EXAMPLE B:
    arcpyproduction.mapping.ClipDataFrameToGeometry (df, ExtentFCLayer, [0])

    del ExtentFC, ExtentFCLayer, ExtentLyr


texts:
Not the first problem I've had with DEMs in QGIS 3.4. Trying to merge 6 DEM tiles together into one. Process for this is simple: go to Raster --> Miscellaneous --> Merge. Once in the menu, select which input files and then select a destination output folder. I get an error that reads like this:

Processing algorithm…
Algorithm 'Merge' starting…
Input parameters:
{ 'DATA_TYPE' : 1, 'INPUT' : ['/Users/nickreale/Documents/Micro Atlas/Chiang Mai/Mapping Data:Info/NorthernThailand/ThaiDEM300_0000240000_0012000000.asc','/Users/nickreale/Documents/Micro Atlas/Chiang Mai/Mapping Data:Info/NorthernThailand/ThaiDEM300_0000480000_0012240000.asc'], 'NODATA_INPUT' : None, 'NODATA_OUTPUT' : None, 'OPTIONS' : '', 'OUTPUT' : '/Users/nickreale/Documents/Micro Atlas/Chiang Mai/Mapping Data:Info/NorthernThailand/northai_1.tif', 'PCT' : False, 'SEPARATE' : False }

GDAL command:
gdal_merge.py -ot Int16 -of GTiff -o "/Users/nickreale/Documents/Micro Atlas/Chiang Mai/Mapping Data:Info/NorthernThailand/northai_1.tif" --optfile /var/folders/42/ld6ywf2d5dsb3g3khrd5j1840000gn/T/processing_69556a0719924167983e1a32cdfe7121/mergeInputFiles.txt
GDAL command output:
/bin/sh: gdal_merge.py: command not found

Execution completed in 0.08 seconds
Results:
{'OUTPUT': <QgsProcessingOutputLayerDefinition {'sink':/Users/nickreale/Documents/Micro Atlas/Chiang Mai/Mapping Data:Info/NorthernThailand/northai_1.tif, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
The following layers were not correctly generated.

/Users/nickreale/Documents/Micro Atlas/Chiang Mai/Mapping Data:Info/NorthernThailand/northai_1.tif

You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.

I've alternatively tried to do this with just two of the 6 DEMs from my folder, thinking that maybe 6 files together were too large.

How do I proceed?


texts: Controlling label/attributes of grid/fishnet index using ArcGIS for Desktop and ArcPy?

texts:
I have used the Grid index feature tool to create an index grid. The label/attribution of the page names are populated in a "row" orientated manner. eg lef to right a1,a2,a3,a4,a5,a6
b1,b2,b3,b4,b5,b6
c1,c2,c3,c4,c5,c6

How can I calculate a field so I can have the page names label/attribution in a "column" orientation?eg:
a1,b1,c1,d1,e1
a2,b2,c2,d2,e2
a3,b3,c3,d3,e3
a4,b4,c4,d4,e4


texts: Editing two layer at same time using QGIS

texts:
How can  Iedit two layer at the same time?

I have node, link layer.

If I edit a node's location, then that should change the link layer automatically.


texts:
I need to preload images (.png) / tiles in my own

source: new ol.source.XYZ

,which can it run fast to move and zoom in mobile devices.

I used ...libs/ol3/3.7.0/ol-debug.js and as urls: urlArray ..

var urlArray = ["/{z}/{x}/{y}.png"];

and tried following parameters to force pre-loading more images than default zoom of map :
'buffer','preload'; but nothink helped/worked. As first I minimize the size of every tile, but it helped not too much.

Can be helped to preload images/tiles throught CSS and

div.none {display:none;background-image:url()}

or exist any JS function to preload images from my folders /{z}/{x}/{y}.png ?


texts: How to export image/imagecollection (raster) as table in GEE?

texts:
I have calculated LST (land surface temperature) of a city, so I have an image. I can export it to my drive as .tiff, but what I need is a table that contain LST value for each longitude and latitude in that city.
I don't know how to export that image to a table (or something like that which contain value of LST in specific long/lat).
How to do this task using Javascript GEE?


texts: Identifying State Parks in the PAD-US dataset (Protected Areas Database)

texts:
I'm looking to map "State Parks," and I found this phenomenal, but complex database. Has anyone ever worked with the PAD datasets before? (or know any easier way to map State Parks?)

In the PAD download there are several included shapefiles, including: fee, easement, marine, proclamation, and designation. It sounds like there is overlap in the datasets. So far my best guess is to use the fee shapefile, filtering on access='open access' and designation='state historic or cultural areas' or 'state park' or 'state recreation area'. But I'm just not 100% sure if I've gotten everything this way. The total number of features returned is 8,065, which is close to the number noted at stateparks.org (8,565 for FY 2017), but it's still 500 features short. When I add in designation='state other or unknown' it brings the total up to 8,447. Although I suppose this number could have (and likely did) go down since the latest release of the data. Still, it would be nice to know the exact query I should be using.

PAD-US Main Page
https://www.usgs.gov/news/mapping-public-lands-united-states

Data Download*
https://www.usgs.gov/core-science-systems/science-analytics-and-synthesis/gap/science/pad-us-data-download?qt-science_center_objects=0#qt-science_center_objects

*I downloaded the "national shapefile"

I'm aware of all the following helpful documentation already:

Data Structure and Attributes
http://www.protectedlands.net/pad-us-data-structure-attributes/

Metadata XML
https://data.usgs.gov/metadata/Science_Analytics_and_Synthesis/5b030c7ae4b0da30c1c1d6de.xml

PAD-US Data Manual:
https://www.usgs.gov/core-science-systems/science-analytics-and-synthesis/gap/pad-us-data-manual

Disclaimer:
“Any use of PAD-US data for showing recreational lands must be careful to show only those lands that have fee ownership and designations that are appropriate to recreational use – the PAD-US ‘Access’ attribute is a key tool for these applications, defining what is Open, Restricted/Permitted, and Closed for public access.”

Citation
U.S. Geological Survey (USGS) Gap Analysis Project (GAP), 2018, Protected Areas Database of the United States (PAD-US): U.S. Geological Survey data release, https://doi.org/10.5066/P955KPLE.


texts: QGIS 2.16: processing.runalg fails when run outside of QGIS in a custom application

texts:
I am developing the RQGIS package. To run QGIS from within R, I call QGIS via the command line. This worked perfectly with QGIS 2.14. Running the same code with QGIS 2.16, however, produces an error.
Here is what I do:

First, I set all necessary paths in the command prompt:

@echo off
SET OSGEO4W_ROOT=D:\osgeo4w_qgis16
call "%OSGEO4W_ROOT%"\bin\o4w_env.bat
@echo off
path %PATH%;%OSGEO4W_ROOT%\apps\qgis\bin
set PYTHONPATH=%PYTHONPATH%;%OSGEO4W_ROOT%\apps\qgis\python;
set QGIS_PREFIX_PATH=%OSGEO4W_ROOT%\apps\qgis
rem open python
python.exe

Subsequently, I run following lines in Python:

import os
import sys
from qgis.core import *
import qgis.utils
from osgeo import ogr
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.gui import *
QgsApplication.setPrefixPath('D:\osgeo4w_qgis16\apps\qgis', True)
app = QgsApplication([], True)
QgsApplication.initQgis()
sys.path.append(r'D:\osgeo4w_qgis16\apps\qgis\python\plugins')
from processing.core.Processing import Processing
Processing.initialize()
import processing
processing.alglist()
processing.alghelp("grass7:v.voronoi")

This works. However, calling processing.runalg...

from processing.tests.TestData import points
result = processing.runalg('grass7:v.voronoi', points(), False, False, '270778.60198,270855.745301,4458921.97814,4458983.8488', -1, 0.0001, 0, None,)

...produces this error message:

Error in sys.excepthook:
Traceback (most recent call last):
  File "D:\osgeo4w_qgis16\apps\qgis\python\qgis\utils.py", line 196, in qgis_excepthook
    showException(type, value, tb, None, messagebar=True)
  File "D:\osgeo4w_qgis16\apps\qgis\python\qgis\utils.py", line 107, in showException
    open_stack_dialog(type, value, tb, msg)
  File "D:\osgeo4w_qgis16\apps\qgis\python\qgis\utils.py", line 142, in open_stack_dialog
    iface.messageBar().popWidget()
AttributeError: 'NoneType' object has no attribute 'messageBar'

Original exception was:
Traceback (most recent call last):
  File "<stdin>", line 10, in <module>
  File "D:\osgeo4w_qgis16\apps\qgis\python\plugins\processing\tools\general.py", line 75, in runalg
    alg = Processing.runAlgorithm(algOrName, None, *args, **kwargs)
  File "D:\osgeo4w_qgis16\apps\qgis\python\plugins\processing\core\Processing.py", line 304, in runAlgorithm
    ret = runalg(alg, progress)
  File "D:\osgeo4w_qgis16\apps\qgis\python\plugins\processing\gui\AlgorithmExecutor.py", line 52, in runalg
    progress.error(e.msg)
AttributeError: 'NoneType' object has no attribute 'error'


texts: Huge offset in Lambert Conformal Conic data using ArcMap

texts:
I want to display my lat/lon climate data, which is in Lambert Conformal Conic. But when I load it in ArcMap (also in QGIS) it will always be displayed with a huge offset. The data is located in northern Germany but is always displayed somewhere in the ocean west of Africa.
The projection has the following parameters. (This is the griddescription generated by my climate data operator - should be accurate).
xunits    = "m"
yunits    = "m"
false_easting = 4000000.f
false_northing = 2800000.f
inverse_flattening = 298.2572f
latitude_of_projection_origin = 52.f
longitude_of_central_meridian = 10.f
scale_factor_at_projection_origin = 0.01745329f
semi_major_axis = 6378137.f
semi_minor_axis = 6356752.f
standard_parallel = 35.f 65.f

I have tried all I can – the data will still show up in the middle of the ocean. I tried transforming and projecting it a million times but that’s of no use. It always comes out the same.
Here is a Link to the actual data set: [Climate Data], https://we.tl/t-M1hidT4Wg4
Does anyone know what’s the mistake I am not seeing here?


texts:
I have two separate computers, each of which are running ArcGIS Desktop 10.3.1 and PostgreSQL and PostGIS 2.1.7. The only difference is one is running PostgreSQL 9.3.4, the other is running 9.3.5.

I am trying to make a simple (or what should be simple) connection to DBMS. I have properly copied all of the PostgreSQL client files ( libeay32.dll, libiconv-2.dll, libintl-8.dll, libpq.dll, and ssleay32.dll) to the ArcGIS bin directory. The pg_hba.conf for each is set to accept client connections. And, of course (of course!!!), I have no problem whatsoever connecting to either using other clients, such as QGIS. However, when I try to add a database connection using ArcGIS Desktop, it works splendidly on one machine (the one running 9.3.4), but not on the other. The error I receive is “Unable to connect to database server to retrieve database list; please verify your server name, user name and password info and try again. Cannot connect to database because the database client software failed to load. Be sure the database client software is installed and correctly configured.”

I have searched high and low to solve this issue and, as a result, have wasted more time than I would like to admit. I am debating if it is more maddening than frustrating or vice versa. If the option was to use another client, such as QGIS, I would. Obviously, it is not. Since I can connect just fine on one machine, I know it actually should work (and the version difference of PostgreSQL really shouldn’t matter).


texts: How to upload a shapefile using WFS service?

texts:
I'm trying to make an application that let the user upload his shapefiles and manage him (show/hide, delete, edit, etc...)  I read in the GeoServer documentation that a WFS is the service that can help me to accomplish my purpose.

I searched enough to find an example on how to upload shapefile using this service but I found nothing so I'll be glad I someone help me out.


texts: Python stops working on Windows when calling layer.GetFeatureCount() of OGR

texts:
My question concerns OGR in Python3 on Windows. So I've installed gdal from the wheel file, and I started playing with OGR functions, but there's a problem I'm encountering that I can't understand. If GetFeatureCount() is called on a layer that was opened in a separate file Python stops working, but not if all my code is in the same script (reading & get count of features).

The code below works:

from osgeo import ogr

filename = '<path>'
driver = ogr.GetDriverByName("ESRI Shapefile")
datasource = driver.Open(filename)
layer = datasource.GetLayer()
print('Features count:', layer.GetFeatureCount())

But if I move the opening of the shapefile's layer to a separate file, it doesn't work anymore (Python stops working somehow):

class IO:
    @staticmethod
    def read_vector(filename):
        # load vector
        driver = ogr.GetDriverByName("ESRI Shapefile")
        datasource = driver.Open(filename)
        layer = datasource.GetLayer()

        return layer

And it called from this file:

from gdal_io import IO

filename =  '<path>'
layer = IO.read_vector(filename)
print('Features count:', layer.GetFeatureCount())

Is something necessary after reading the shapefile vector with OGR?


texts:
I counted points in polygons in PostGIS like this

UPDATE hexagons SET count = (SELECT count(*) FROM points WHERE ST_Within(points.geom, hexagons.geom));

That came back with that

UPDATE 1554

That worked fine. It updated the new created column with the amount of points inside each polygon (6.6 mio points, 50k polygons).
I started up QGIS to display the distribution of my points.

right click on the layer
properties
style -> graduated -> picked my colorschema and so forth

finally, I clicked okay and it looked like that (underlaying my point layer with which I did the 'Count Points in Polygon' in PostGIS with

Well, I thought maybe something went wrong with my querry in PostGIS so I checked how many rows there are empty in my polygon layer. It turns out just two rows are empty.

What is weird though is that in the scale it only goes until 999 while there are rows which go well over 100000. I think that might be the problem.
Whats going wrong here? Does QGis have problems with rows where the value is > 1000?
I cannot find anything to that.


texts: How to free the memory of return value of OGRGeometry::Buffer()?

texts:
I have in C++ (gdal 2.3.2)

OGRPolygon poly;
OGRGeometry* geom = poly.Buffer(1e-12, 1);

now how do I free the memory allocated for geom after done using it?


texts: Spatial+Temporal database request for 4D intersection

texts:
Here is the function I need to implement. My application deals with mobile objects with variable position and range (think of radar range). As I can't test each minute if n mobiles detect each other (server would be overloaded), I need to precalculate "future" position and time when they'll detect each other from their planned route and only register the "future" event in a stack. Each time a user changes the planned route, I'll recalculate event related to all intersecting routes from other mobiles.

The best way to do that is to make a kind of 4D intersection request and update, and be able to detect when and where mobiles will "collide".

Is it possible with a classic geodatabase like PostgresQL+PostGis and another plugin, or do I need to look at specific n-dimensional database ?


texts: Creating Latitude & Longitude from polygon shapefile using QGIS

texts:
Could you help me to create long & lat (x and y coordinates in degrees) for polygon shapefile?

When I used directly QGIS export vectors to csv file, it didn't work.

Here is the example file: https://drive.google.com/drive/folders/1lfkwMQADzsiMKCZsOeuFqQAHGUE9Qni8?usp=sharing


texts: how to remove and change qgsvertexmarker in QGIS using c++ code?

texts:
In QGIS following code is used to set marker at point A.

  QgsVertexMarker * marker = new QgsVertexMarker(mPlugin->iface()->mapCanvas());
  marker->setColor(Qt::black);
  marker->setIconType(2);
  marker->setCenter(A);

How to delete this marker using C++ code?

In python it is

 canvas.scene().removeItem(marker)

I am using this marker in a loop where point A will change continuously so i would like to delete marker at previous position and insert in new.

when i use below code marker is deleted permanently.

  mPlugin->iface()->mapCanvas()->scene()->removeItem(marker);

so how to use this marker to plot moving point?


texts: Changing default parameters of "Merge Selected Entities" in QGIS 2.18

texts:
I am using QGIS 2.18, and I am looking for a way to change default values of the Merge Selected Entities algorithm.
To make it short, I have in my shape file buildings that are divided in several small parts (whereas it should not), and I would like to merge those part to get back a single building. However, one of my attributes is the height, and the default value of the merging algorithm is to pick one height among the merged entity, whereas I would like it to be always the maximum (I have quite a lot of buildings to fuse and I would like to avoid changing parameters every time I need to merge entities). Is there a way to change this through some files?


texts: How to customize the visibility of buttons in QGIS?

texts:
The built-in toolbar named "Composer" has an icon "Duplicate Composer".

I wish to toggle the visibility of this button in the toolbar.

How do I achieve this using pyqgis?


texts: How can I style many layers of vector grid once in QGIS?

texts:
After I created many of vector grid layers, I now need to style them in a batch system. Is there any way to do this?

If that is not possible, how can I use the layer duplicate option to inherit the style from the original grid; and then change the extents of the duplicated layer?


texts: How to resize the symbol when the map extent gets changed or any event is triggered

texts:
I need to resize the PictureMarkerSymbol when the map extent gets changed.
I am able to catch the event with watch() but somehow resizing does not apply to the map.
I've tried below codes but it does not work.

<!DOCTYPE html>

<meta charset="utf-8">
<meta name="viewport" content="initial-scale=1, maximum-scale=1, user-scalable=no">

<link rel="stylesheet" href="https://js.arcgis.com/4.2/esri/css/main.css">

If I call resizeBuildingImg() function directly from view.then(.... ), not from watch function, it works as expected. So I think I miss something. But not sure what I don't get on this api.


texts: Understanding backslash in arcpy.SelectLayerByAttribute syntax?

texts:
I copied this line from Python code, but I don't understand the \" \" around the word District.

arcpy.SelectLayerByAttribute_management(Community_areaFC,
     "NEW_SELECTION",
     "\"District\" = 3")


texts: Virtual Raster Mosaic - gdalbuildvrt - data gaps

texts:
I'm trying to create a virtual mosaic of several raster images (.tif) in R using the function "gdalbuildvrt" from the "gdal" package. The raster images are flight lines from a drone flight over a field.

I specified the corresponding NoData value for the input raster in the function

gdalbuildvrt(gdalfile = paste(--path--),
         output.vrt= paste(--path--),
         separate=FALSE, srcnodata="-3.4e+38", vrtnodata="-3.4e+38",
         overwrite = TRUE, r="average")

The "NoData" areas of the individual raster files seem to overwrite each other, resulting in a mosaic with missing data.

The raster images are in formatted "Float32 - Thirty two bit floating point" data type and the same value specified as their NoData value (I checked in QGIS, checking the value from an area without data).

Manually merging the raster files in QGIS using Raster --> Tools --> Merge is successful and yields a correct mosaic without the white overlapping areas.

gdal_merge.bat -ot Float32 -of GTiff -o
C:/.../OUTPUT.tif --optfile
C:/...\mergeInputFiles.txt

Also, when changing the data type to 8bit (0 to 255) beforehand and defining "255" as srcnodata value, the resulting mosaic from the gdalbuildvrt function is correct. However, that's not an option since that requires stretching the data range of the individual rasters - which I don't want to do before the whole mosaic is created.

What is the fault with my approach? I've also tried specifying various other nodata values for srcnodata/vrtnodata ("nan", "0", "3.4e+38") but I get the same results.


texts: Tool for downloading map data from Public websites which has Google as Base map

texts:
I am trying to download the map data (coverage map) from the MAXIS web site. I couldn't find any option to download it directly. This web service has coverage map of two different band widths with Google as base map.

Can we take this URL https://maxisnetwork.tuxgeo.com/v2/) and put it in some applications and download the coverage map data?

I have done some research but couldn't find the solution.


texts: ArcMap Field Calculator gives ERROR 999999

texts:
I am trying to do a fairly simple calculate but I keep getting ERROR 999999. I am running ArcMap 10.7.

My Input is:

(500,000/16) + ( (500,000/25) * [Join_Count])


texts: MySQL Using ST_CONTAINS to obtain points within a bounding box

texts:
Can someone give me an insight into why the following MySQL code is not working.  I have taken some of this code from this link and it is all about just proving the concept before applying to my site.  What I want to do is find the points on a line that are bounded by a polygon.

The line is defined by:

LINESTRING (-4 5,5 5,10 5,12 5)

and is entered into the mySQL new_table row new_tablecol as a BLOB.   If I inspect the field it is consistent with the LINESTRING above.

I now have defined a bounding box and SQL using the following syntax:

SET @g1 = GeomFromText('Polygon(0 0, 10 0, 10 10, 10 0, 0 0)');
SELECT * FROM `new_table` WHERE ST_CONTAINS(`new_tablecol`,@g1)

Clearly the bounding box includes some of the LINESTRING, but no matter how hard I manipulate this I can't get any points returned.

Can anyone help?


texts: Deriving Topographic Information from Geotagged Photos?

texts:
I wish to extract elevation data from geotagged photos?

I have worked with DEM's and the open source CLI program: GMT for my graduate studies. I own a DJI Phantom 3 and used the trial version of Drone Deploy.

They are able to take geotagged jpg photos and create elevation data from it but how is this possible?

In the trial version they do not allow direct download of raw elevation data, but they do supply a Geotiff of the survey area, and using GMT I can convert the geotiff to netCDF file where then I can plot it again. However the elevation data isn't accurate as it would be on the DroneDeploy website itself.

Therefore, I am forced to work my ways backwards to find how to extract elevation data from geotagged jpgs.


texts: add/show layer(group) on click [leaflet]

texts:
I am a beginner in leaflet.

Basically, I have a circle. I would like to show another circle when clicked on this circle. And hide it again when clicked for the second time. How is that possible?

 //the circle which will haven on click event
    var circle1 = L.circle([44.6663888888889, 6.775],{radius: 150, color:'white',weight:.5, opacity:1,fillColor: '#ff9900',fillOpacity:1});

//the circle to be added when one click on circle 1
    var circle2 = L.circle([44.6502777777778, 6.79333333333333], {radius: 150, color:'white',weight:.5, opacity:1,fillColor: '#008000',fillOpacity:1});


texts: Seeking GIS Program With Dynamically Updating Layers?

texts:
I am looking for a GIS platform which supports dynamically updating layers.  These could be updated by file, socket, plugin, or any other measure, that isn't important.

The goal is to be able to draw shapes on a map and update those shapes via an external program, with the results being sent (or written to file) to a GIS application which would update its layers accordingly.  Instant updates aren't vital, if the GIS application checks on a timer that is acceptable.

I've spent some time searching and poking through the docs on things like QGIS and I haven't found anything definitive that supports this functionality.  The "monitor layer" option on delimited layers in QGIS is along the lines of what I want, I think, but the map will only update upon a manual action from what I can tell.

Dynamic updating of layers through an external process is the most important requirement for this small project and far outweighs other GIS capability considerations.


texts:
I am wondering is there any opensource UML to Postgresql Mysql tables conversion software which also includes the geospatial datatypes ?

Note:I want conversion from UML to database. and NOT visualization of existing database

I know enterprise architect but it is commercial.


texts: QGIS 2.0 - can't select items in Print Composer

texts: I am having an issue where I can create a print Layout in Print Composer, but once I have clicked away from an item (like map, legend, arrow, etc) I cannot select it again.


texts: Why does Pyproj fail converting lon lat alt to ECEF?

texts:
I'm trying to convert a WGS84 lon, lat and alt to ECEF and Pyproj is giving me really weird results. My code is simply:

lon, lat, alt = 49.74761271, -113.2179781, 0
ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum='WGS84')
lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum='WGS84')
x, y, z = pyproj.transform(lla, ecef, lon, lat, alt, radians=False)
print x, y, z

which prints

inf inf 0.0

Does anyone know what is wrong?

If I implement it myself it works:

def lla_to_ecef_alt(lat, lon, alt):
    rad = np.float64(6378137.0)        # Radius of the Earth (in meters)
    f = np.float64(1.0/298.257223563)  # Flattening factor WGS84 Model
    np.cosLat = np.cos(lat)
    np.sinLat = np.sin(lat)
    FF     = (1.0-f)**2
    C      = 1/np.sqrt(np.cosLat**2 + FF * np.sinLat**2)
    S      = C * FF

    x = (rad * C + alt)*np.cosLat * np.cos(lon)
    y = (rad * C + alt)*np.cosLat * np.sin(lon)
    z = (rad * S + alt)*np.sinLat
    return x, y, z

What's the difference between them?


texts: Python GDAL: Save array as raster with projection from other file

texts:
I have an array of data, and for each datapoint I know the latitude and longitude. I'd like to save it as a GTiff with the same projection as other rasters I have. This is what I've tried so far, but no luck.

import numpy as np
import gdal
from gdalconst import *
from osgeo import osr

def GetGeoInfo(FileName):
    SourceDS = gdal.Open(FileName, GA_ReadOnly)
    GeoT = SourceDS.GetGeoTransform()
    Projection = osr.SpatialReference()
    Projection.ImportFromWkt(SourceDS.GetProjectionRef())
return GeoT, Projection

def CreateGeoTiff(Name, Array, driver,
                  xsize, ysize, GeoT, Projection):
    DataType = gdal.GDT_Float32
    NewFileName = Name+'.tif'
    # Set up the dataset
    DataSet = driver.Create( NewFileName, xsize, ysize, 1, DataType )
            # the '1' is for band 1.
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection( Projection.ExportToWkt() )
    # Write the array
    DataSet.GetRasterBand(1).WriteArray( Array )
    return NewFileName

def ReprojectCoords(x, y,src_srs,tgt_srs):
    trans_coords=[]
    transform = osr.CoordinateTransformation( src_srs, tgt_srs)
    x,y,z = transform.TransformPoint(x, y)
    return x, y

# Some Data
Data = np.random.rand(5,6)
Lats = np.array([-5.5, -5.0, -4.5, -4.0, -3.5])
Lons = np.array([135.0, 135.5, 136.0, 136.5, 137.0, 137.5])

# A raster file that exists in the same approximate aregion.
RASTER_FN = 'some_raster.tif'

# Open the raster file and get the projection, that's the
# projection I'd like my new raster to have, it's 'projected',
# i.e. x, y values are numbers of pixels.
GeoT, TargetProjection, DataType = GetGeoInfo(RASTER_FN)
# Meanwhile my raster is currently in geographic coordinates.
SourceProjection = TargetProjection.CloneGeogCS()

# Get the corner coordinates of my array
LatSize, LonSize = len(Lats), len(Lons)
LatLow, LatHigh = Lats[0], Lats[-1]
LonLow, LonHigh = Lons[0], Lons[-1]
# Reproject the corner coordinates from geographic
# to projected...
TopLeft = ReprojectCoords(LonLow, LatHigh, SourceProjection, TargetProjection)
BottomLeft = ReprojectCoords(LonLow, LatLow, SourceProjection, TargetProjection)
TopRight = ReprojectCoords(LonHigh, LatHigh, SourceProjection, TargetProjection)
# And define my Geotransform
GeoTNew = [TopLeft[0],  (TopLeft[0]-TopRight[0])/(LonSize-1), 0,
           TopLeft[1], 0, (TopLeft[1]-BottomLeft[1])/(LatSize-1)]

# I want a GTiff
driver = gdal.GetDriverByName('GTiff')
# Create the new file...
NewFileName = CreateGeoTiff('Output', Data, driver, LatSize, LonSize, GeoTNew, TargetProjection)

But this results in the following error message:

File "TES2GtifBBB.py", line 25, in CreateGeoTiff
  DataSet.GetRasterBand(1).WriteArray( Array )
File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal.py", line 1082, in WriteArray
  return gdalnumeric.BandWriteArray( self, array, xoff, yoff )
File "/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site packages/osgeo/gdal_array.py", line 256, in BandWriteArray
  raise ValueError("array larger than output file, or offset off edge")
ValueError: array larger than output file, or offset off edge


texts: How to remove ArcGISTiledMapServiceLayer completely

texts:
I have a problem related to "ArcGISTiledMapServiceLayer".
I have an ArcGISTiledMapServiceLayer map service, which I have added to Map object. This service work as BaseMap. Now I want to add another ArcGISTiledMapServiceLayer map service layer to the same Map object. Before adding the new layer, I just removed the previous tiled map service layer by using method objMap.removeAllLayers(). But this method only removes the layer not the extent and it create problem when I add a new layer.

So I want to remove completely the ArcGISTiledMapServiceLayer (layer, extent, scale etc.).

I have tried all the below methods

objMap.removeAllLayers();
objMap.removeLayer(objMap.getLayer('myLayerID'));

Any solution?


texts: QGIS: OGR error when editing Geopackage layers : database is locked

texts:
I started to draw new polygons in a layer within an existing GeoPackage that I previously created.
However, after editing several layers, I sometimes (apart that it doesn't happen when editing one single layer, it looks like a random phenomena) cannot save the layer and I get the error:
Errors: ERROR: feature(s) not added.

Provider errors:
  OGR error creating feature: failed to execute insert : database is locked

The polygons drawn are then lost, which is pretty frustrating.
I googled a lot and found this rather old fix which is apparently already applied in QGIS 3.12 (my version).
Any suggestion?


texts: Getting Error 999999 from arcpy.SpatialReference.name?

texts:
I'm trying to build a set that consists of spatial reference names using .name method from the arcpy.SpatialRefernce() class. When I do this for a single prj file there is no issue, but when I try to do it for a set of prjs as a loop it returns a 999999 error:

import arcpy

#overwrite existing output from previuous trials
arcpy.env.overwriteOutput = True

#create output file location through arcpy.env.workspace
arcpy.env.workspace = r"C:\RS_Data\Python\Geog485\Lesson2"

#designate shapefile that has the desired projection
outPRJ = r"C:\RS_Data\Python\Geog485\Lesson2\StateRoutes.shp"

#call on the prj file to get the spatial reference object and name
targetPrj = outPRJ[:-3]+"prj"
SpRef = arcpy.SpatialReference(targetPrj)
SpRefName = SpRef.name
print SpRef.name

#create list of fc prjs and remove .shp file names from previous trials
x = "_projected.shp"
inFiles = arcpy.ListFeatureClasses()
for file in inFiles:
    if x in file:
        inFiles.remove(file)

inPuts = set()
for file in inFiles:
    file = file[:-3]+"prj"
    inPuts.add(file)
    print file

for prj in inFiles:
    sr = arcpy.SpatialReference(prj)
    if sr.name != SpRef:
        inPuts.add(sr)

here is the error:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec cmd in globals, locals
  File "N:\Python\Reproject.py", line 1, in <module>
    import arcpy
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\arcobjects\mixins.py", line 927, in __init__
    self._arc_object.createFromFile(item)
RuntimeError: ERROR 999999: Error executing function.


texts: Are there any open source tools which can be used for map composing, preferably GUI based. I use Quantum GIS sparingly but I am interested to know if there are others. The tool must be able to read raster/vector layers in standard formats and allow exporting in various formats. Something like GeoPDF export would be a great plus.


texts: What kind of input rasters should I use in a Spatial Analyst Corridor process?

texts: I´m trying to map corridors between stands of hardwood. I´ve tried cost path, almost all the time get strange or empty outputs. Now I saw the other tool called Corridor and thought I might give it a shot. But what should I have as inputs? It makes sense to me to use a cost distance raster as one of the inputs but the other? And what does a good output looks like - (so I know if this is the right tool for me)? Tried with the ArcGIS help and online guides but they are very unspecific if you don´t want to go into programming.


texts: Representation of road networks at junctions when line width is varying with attribute

texts:
I have a road network and for every road segment the average number of daily traffic. I need to draw a map where the line width for each road segment is proportional to the daily traffic. I can do this in ArcMap using the proportional symbols symbology, or in QGIS using a data defined property in the line style (where the pen width is proportional to the daily traffic using some scaling formula). However, at junctions, the graphic representation is not optimal because the line ends are (of course) exactly overlapping (like in A in the image below). A more correct representation would be that the lines are adjacent at the junctions, like in B in the image below.

I could manage some part of this problem using the date defined properties in QGIS on the line offset, but it is still far from being optimal. It also gives the new problem of misalignment at the end of the segment (top right in the image below), and a problem with the line labels. This is for now the best I could manage:

Is there any way to achieve this type of rendering like in B in the first image, in QGIS or ArcGIS/ArcMap? Or is there any other GIS software which gets this type of rendering correct?


texts: Extracting values from raster to get elevation values in R?

texts:
I am trying to generate an isoscape using IsoriX in R. Here how my data looks:

> head(dataC)

  siteID  long   lat       animalID       species tissue.value

1      A -5.88 53.82  CC1_33.2/1_15             C       -18.46

2      A -5.88 53.82  CL1_33.2/1_15             L       -19.98

3      A -5.88 53.82  CT1_33.2/1_15             P       -19.17

4      A -5.88 53.82  STA_33.2/1_15             S       -18.49

5      A -5.88 53.82 MESO_33.2/1_15             M       -27.05

6      B -6.03 53.41    CC1_43/3_15             C       -17.11

> str(dataC)

'data.frame':   322 obs. of  6 variables:

 $ siteID      : Factor w/ 38 levels "A","AA","B","BB",..: 1 1 1 1 1 3 3 3 5 5 ...

 $ long        : num  -5.88 -5.88 -5.88 -5.88 -5.88 -6.03 -6.03 -6.03 -5.87 -5.87 ...

 $ lat         : num  53.8 53.8 53.8 53.8 53.8 ...

 $ animalID    : Factor w/ 322 levels "AA1_36/8_15",..: 138 187 257 318 267 144 272 294 189 273 ...

 $ species     : Factor w/ 7 levels "A","C",..: 2 3 6 7 4 2 4 5 3 4 ...

 $ tissue.value: num  -18.5 -20 -19.2 -18.5 -27.1 ...

However, as the elevation is unknown at the sampling sites, it says it has to be extracted from a high resolution elevation raster using the function extract. However I cannot make it work...

The area I need to get on is the Irish Sea:

>map("worldHires",xlim=c(-6.7,-2.7), ylim=c(52,56), col="gray90", fill=TRUE,
myborder = 0.0001)

This is what I am typing:

## Creating Raster

>library(sp)  # classes for spatial data
>library(raster)  # grids, rasters
>library(lattice)
>library(latticeExtra)
>library(RColorBrewer)
>library(rasterVis)  # raster visualisation
>library(maptools)
>library(rgeos)
>library(dismo)
>library(XML)
>library(rgdal)

>mymap <- gmap("England", type = "satellite", filename = "England.gmap")
>plot(mymap)

>select.area <- drawExtent()
### now click 2 times on the map to select Irish Sea
>mymap <- gmap(select.area)
>plot(mymap)

### Save the map as a file in your working directory for future use
>mymap <- gmap(select.area, filename = "IrishSea.gmap")

## Get elevation values

### Read point data, and convert them into spatial points data frame
>pointCoordinates <- read.csv(file.choose(),h=T)
>coordinates(pointCoordinates) <- ~ x + y

### Extract raster value by points
>rasValue <- extract(mymap, coordinates(pointCoordinates))

### Combine raster values with point and save as a CSV file
>combinePointValue <- cbind(pointCoordinates,rasValue)
>write.table(combinePointValue,file= "combinedPointValue.csv",
            append=FALSE, sep= ",", row.names = FALSE, col.names=TRUE)

--


texts: I am trying to create a path from the blue square to all the green dots. The blue square represents an existing road and the green dots represent centroids of forest stands. The total path distance shall be as short as possible and the path shall be the least cost-path (based on Slope Raster). The least cost path can be done with ArcGIS Cost Path. But I don't know how I could create the shortest possible path.


texts:
I would like to find the minimum, maximum and (weighted?) average value of the raster cells intersected by a line (vector) layer. This is similar to what common "Zonal Statistics" tools do, but they all work on polygons only, not lines.

It appears that Hawths Tools have something similar called Line Raster Statistics, but I'm looking for an Open Source (non esri) solution.

I use QGIS, so either a native solution, plugin, or SQLite or PostGIS solution would be preferred.


texts: Assigning z values to polygon feature vertices in ArcMap

texts: I am trying to assign z values to vertices of rectangular polygon features in ArcGIS 10.5 by using a DEM raster. I do not want to assign min, max, or mean z values to the polygons themselves, but I am eventually trying to visualize my polygon as a TIN-like surface in another program (EVS) without producing a TIN; the final result must be a polygon. I've tried the Interpolate Shape tool, but it seemed to assign z values to points within the rectangular polygon features as well, creating a very large and unmanageable shape.


texts: Working with Google Maps/ArcMap lat/long

texts:
I know a bit about lat/long, projected and geographic coordinate systems etc., and a bit about geodesy and datums.

I have a lat-long in deg/min/sec provided to me by someone else (point location = northern Italy). I have a map in ArcGIS in the same region and I'm constructing a grid over it for some analysis, centred on the lat-long point above.

In ArcGIS, I'm working with projection Roma 1940 Gauss Boaga Ovest - the data frame and map layer are both in it - and when adding the XY point also in this projection, it ends up pretty much where it should be (but not exactly)

The thing is, when putting these coordinates into Google Earth, the point is around 25 or 30 metres south east of where it is in ArcGIS. The provider is saying the Google Earth location is where it should be.

So, where does the difference derive from?

I understand Google Earth is a WGS84 geographic coordinate system so I tried a point file in WGS84 and added the XY and it still ends up in the slightly incorrect position in ArcGIS as outlined above. Also have tried various combinations of WGS84 and Roma 1940 Gauss Boaga Ovest and can't get it right.

Would the Roma 1940 datum skew it slightly or is it nothing to do with that?

I see that Google Maps actually uses WGS84 Web Mercator projected coordinate system, which has its own problems.

Would ArcGIS actually be more accurate then?

I asked the provider and they said the coordinates provided have come from Google Maps.

I've tried to reproject etc. and nothing will change the position from the red dot (see left of image). I've since discovered Google Maps is a WGS_1984_Web_Mercator projection, based on a sphere. I've attempted to create a custom transformation of datum in ArcMap from 'GCS_WGS_1984_Major_Auxiliary_Sphere' to 'GCS_Roma_1940', no change to dot position when entering the coordinates.

As you can see, the Google Maps on the right is in a slightly different position, by roughly 35m, to the ArcMap image on the left. I'm pretty much out of ideas.


texts: Latlng to UTM conversion with forced zone

texts:
I'm trying to convert a set of latlng coords to UTM but the points in the set are across two different UTM zones, 32T and 33T.

This causes two near points like

44.51601173544275, 11.992860160820502
44.51601173544275, 12.012860160820502

to be converted to

32T 737655.84708598 4933706.271726263
33T 262397.6792462121 4933704.362693829

Code is:

LatLng latlng = new LatLng(44.51601173544275, 12.012860160820502);
latlng.toWGS84();
UTMRef utm = latlng.toUTMRef();
System.out.println(utm);
// Prints 33T 262397.6792462121 4933704.362693829

For some reason, when qGIS reads those points seems to ignore the zone and just read the numerical coords, so the final result is that the points are very distant from each other:

Black lines are 32T and 33T.

So, the question is: is this a qGIS bug or a conversion issue? Can I force the conversion to use a specific zone (32T) as a workaround? Right now I'm using Jcoords for conversion but in my project I already have Geotools loaded.

Thanks, regards
Maurizio


texts: Can 'Source' and 'Target' in the roads table for routing be duplicate?

texts:
I am using pgr_ksp for finding multiple paths for single source and destination. I am wondering if 'source' and 'target' in a roads table are unique or not. because when i am doing routing on my road network table and displaying their geometry on the basis of their 'target' value, it is giving me extra link roads everywhere. This is the query i am using

   SELECT pt.gid, seq AS route, id1 AS node,  id2 AS edge, cost, geom as line
   FROM pgr_ksp(
   'SELECT gid as id, source::int4, target::int4, length::float8 as cost FROM my_table',
    2087, 1428, 2, false
    )
    as di JOIN my_table pt
    ON di.id2 = pt.target

The red color is showing the output of the query.
I know there is nothing wrong with the routing result but when i am trying to get the geometry, 'target' value is i think the only way to get it for display

I want the result as i am showing in the 2nd image with Red color, but right now i am getting extra road segments (shown in Green color) because i am choosing to display geometry from target value.

Kindly suggest me what to do in order to get better result


texts:
Considering that depth denotes the distance below sea-level to a point on an object located below sea level, and considering that height denotes the distance to some point above the surface of the earth as measured from a point on surface of the earth, it seems there ought to be a term that denotes the distance to some point below the surface of the earth as measured from a point on the surface of the earth.

For example, suppose someone dug a mine shaft down into the earth from an elevated point, but not so far down that it descends below sea-level; what term ought we use to denote the points below the surface of the earth?


texts: Use Hibernate with Oracle and Esri St_geometry type

texts:
Can anyone tell me whether I could use Hibernate with esri ST_GEOMETRY type?
I have found that Hibernate Spatial does support SDO_Geometry, but nothing for ST_Geometry.

I am using oracle 11g and arcgis 10.2


texts: Connect points along bathymetry contours

texts:
I have a table with start and end point data in each row and need to convert them into lines, but I would like something other than a straight segment.
Ideally, I want to connect the points via a segment that is determined by the bathymetry of the area (so along the same bathymetry as the start and end points). I have both raster bathymetry layers and contour shapefiles.
The only way I can think to do this is create a raster layer with the difference in bathymetry from the mean bathymetry of start and end points and create a least-cost path with large differences being very costly.
I wanted to float the idea out there though to see if there are any automated or different ways to do this.
Here is a diagram of what I have/want:

The black lines are bathymetry contours, and the red points are examples of my start and end points. The blue line is the straight line segment that is how ArcMap tools would connect the two points.
What I need is a line connecting them that follows the bathymetry contours, i.e. the path with the least amount of change in bathymetry from the start and end points.

After exploring the idea, it is not feasible to simply snap the points to the bathymetry contour shapefile and split the lines, that would change the location of the points too much.


texts: How can I detect if my linestring is near an OSM node?

texts:
If I have a linestring of objects (different from the road OSM XML file) and an OSM XML file that contains nodes and ways (extracted down to roads only), how can I detect if my linestring passes through a node?

This is given that I have already made a script that has the list of intersecting nodes from ways, so I just basically want to bring it down to the most significant ones, where my linestring passes through these intersecting nodes.


texts: QGIS 3 does not show geometry tables with NULL geometries

texts:
I just noticed that in QGIS 3 (3.0.1, Lubuntu 16.04 LTS), QGIS does not show tables in which the geometry column has all NULL values.

For example, I took a table tt with a geometry column geom, and

UPDATE tt SET geom = NULL;

Then if I refresh the QGIS connection, the table disappears. The 'geom' column looks as follows in \d+ tt:

                                     Table "public.tt"
 Column |   Type   | Collation | Nullable | Default | Storage  | Stats target | Description
--------+----------+-----------+----------+---------+----------+--------------+-------------
...
 geom   | geometry |           |          |         | main     |              |

Has anyone had similar issues?


texts: Using pyQGIS's iface or QgisInterface modules in a standalone script?

texts:
According to this thread on using the iface module in a standalone pyQGIS script it is not possible to use the iface module in a standalone script, only being accessible within a QGIS Plugin.

In the official QGIS documentation for the QgisInterface module it states that it is accessible only via QGIS Plugins and not standalone scripts

My question is how I could access the map canvas and the manually set extent of a QGIS project using the pyQGIS module or set the active layers in a given QGIS project as most of the examples use the qgis.utils.iface.setactivelayer() function hasn't worked for me up until now


texts: How to use third-party package in a QGIS plugin?

texts:
I'm developing a QGIS plugin and I would like to use arrow or requests packages as I'm used to manipulate dates and REST APIs.

QGIS does not come with these packages (at less arrow), so when I try to access my plugin in QGIS I've got :

  "ImportError: No module named arrow"

How can I specify my requirements?


texts: How can I get GeoServer WMS layer attributes in QGIS?

texts: I am trying to connect GeoServer wms layers in QGIS for the purpose of web mapping. The layers are shown properly in QGIS but the layers attribute that means layer information is not visible in layer properties.


texts: How to get neighborhood values count using a geoprocessing tool in ArcGIS?

texts: I have a raster (black=0, red=1, and white=NoData). I am wondering if there is any tool in arcgis which can help me only count number of 0s and 1s separately in a 3*3 neighborhood. I need an attribute table with that information. Any hint of relevant tool would be appreciative.


texts:
I'm using Google Earth Engine with the Python API. I have this piece of code:

centroids = imNDVI4.reduceToVectors(
        geometry = test2,
        geometryType = 'centroid',
        scale = 30,
        crs = 'EPSG:4326',
        eightConnected = True
    )

I want to save the centroids within this variable, into a TXT file as a list of coordinates.
How can I do this?


texts: Specific Spatial Analyst tools failing while running extension using ArcGIS server

texts:
I am having having trouble running ExtractValuesToPoints on a remote server machine running ArcGIS Server (with ArcInfo and Spatial Analyst installed). It fails and throws a standard COM error. I have checked for any obvious issues with inputs/outputs and it runs just fine on several desktop machines running ArcGIS 10 Desktop (with ArcInfo and Spatial Analyst).

The interesting thing is the Spatial Analyst tool Con works just fine in a previous geoprocessing step, which leads me to believe that it's not a licensing issue. I have even tested the inputs/outputs on the remote machine using the ArcGIS UI and it completes.

I apologize for the vagueness of the question, but I really can't seem to find any reason for one spatial analyst tool to work and another not (other than obvious issues with the data). I'd be happy to fill in any gaps if more information is needed. I guess I'm looking for a list of things I should be checking that I may have glanced over.


texts: ArcGIS 10.1 - Using Polygon Neighbors with dissolved polygons

texts:
I've used polygon neighbors to find all neighbors of a polygon layer, which I then use to look up all possible neighbor combinations and dissolve the polygons into all possible combinations of two polygons. I'd then like to run this again and build all possible three-polygon combinations.

However, when I run the tool again on the dissolved polygons layer it creates an empty table. I know the data is fine - is there another step i need to take, or is there another tool I could use to achieve this? Thanks in advance.

EDIT:

This is the script that i use after running the polygon neighbors tool to create the dissolves and then append them to 'Sites'. I then manually run polygon neighbors again on 'Sites', but get an empty table:

arcpy.env.workspace = "C:\Users\Ant\Documents\ArcGIS\ITN.gdb"
FC = "BLPUs"
TABLE = "BLPUs_PolygonNeighbors"
DISSOLVED = "Sites"
FIELDS = ["src_OBJECTID", "nbr_OBJECTID"]
with arcpy.da.SearchCursor(TABLE, FIELDS) as cursor:
    for row in cursor:
        OID1 = str(int(row[0]))
        OID2 = str(int(row[1]))
        arcpy.MakeFeatureLayer_management(DISSOLVED, "DISS_LYR")
        OUTPUT_NAME = "Dissolved_%s_%s"%(OID1, OID2)
        SQL = "OBJECTID = %s OR OBJECTID = %s"%(OID1, OID2)
        arcpy.SelectLayerByAttribute_management("DISS_LYR", "NEW_SELECTION", SQL)
        arcpy.Dissolve_management("DISS_LYR", "DISS_OUT")
        expression1 = arcpy.GetCount_management("DISS_LYR")
        fieldName1 = "BLPUCount"
        fieldName2 = "OID1"
        fieldName3 = "OID2"
        fieldName4 = "OID3"
        fieldName5 = "OID4"
        fieldName6 = "OID5"
        arcpy.AddField_management("DISS_OUT", fieldName1, "LONG")
        arcpy.AddField_management("DISS_OUT", fieldName2, "LONG")
        arcpy.AddField_management("DISS_OUT", fieldName3, "LONG")
        arcpy.AddField_management("DISS_OUT", fieldName4, "LONG")
        arcpy.AddField_management("DISS_OUT", fieldName5, "LONG")
        arcpy.AddField_management("DISS_OUT", fieldName6, "LONG")
        arcpy.CalculateField_management("DISS_OUT", "BLPUCount", expression1)
        arcpy.CalculateField_management("DISS_OUT", "OID1", OID1)
        arcpy.CalculateField_management("DISS_OUT", "OID2", OID2)
        arcpy.CopyFeatures_management("DISS_OUT", OUTPUT_NAME)
        arcpy.Delete_management("DISS_OUT", "DISS_LYR")

fcList = arcpy.ListFeatureClasses("Dissolved*")
for fc in fcList:
    arcpy.Append_management(fc, DISSOLVED, "", "", "")


texts: ArcSDE how to Change Allow null values to "NO"

texts:
I have created a few feature classes on an SDE using ArcCatalog and initially when I created them I set "Allow null values" to yes. I was wondering if this can be changed back to "No"? It is ArcSDE 10.0 on Oracle 11g. I need to disallow null values to minimise the amount of time spent on QC of null values inserted by default by the user.

See attached screenshot


texts: Packing Polygons within polygon using ArcGIS Desktop?

texts:
I have a Boolean raster.

In the grey areas of the raster I would like to fit a given size polygon within a contiguous extent.

Basically, I have an irregular polygon, and i would like to "fit" a known polygon within the extent of the irregular polygon as many times as possible.

Direction of the polygon does not matter, and it could be a square.  I would like for it to fit graphically, but if it just attached a number to the polygon (# that fit) that would work as well.

I am using ArcGIS Desktop 10.


texts: Extract values from rasterstack TIFF for point data locations in Python

texts:
I am working on machine learning algorithms for predicting wildfire sites on state level. I have a CSV of wildfire instances, that has 37 columns including a geometry column, a lat column and a lon column. I also have a stacked raster, in TIFF format of various climate and other associated variables. I am trying to use Python to extract the values from the stacked raster at the point locations (from the wildfire CSV) and create a data frame.
I have done a similar task before in R, using:
training <- raster::extract(raster_stack, wildfire, df=TRUE)

I have tried a couple of methods, similar to this response Python - Extract raster values at point locations
and this one
with fiona.open('wildfire_gdb/wf_clip.shp') as shp:
    for feature in shp:
        wf_category = feature['properties']['FIRENAME']
        coords = feature['geometry']['coordinates']
        with rasterio.open('/Users/marie/Desktop/GI3/projekt/stacked.tif') as src:
            value = [v for v in src.sample([coords])][0][0]
        raster_values.loc['FIRENAME']=value

However, it just returns a CSV with a single cell of -3.4e+38, and think I am going wrong with calling the geometry or lat/lon values somewhere on a feature, 'FIRENAME' as I am not 100% sure what was meant by 'properties' and 'coordinates' in the code I tried to edit for my own use.
Sorry, I am new to this and not sure what other information you need that would be helpful and have searched everywhere and found nothing so far that gets me the output expected.


texts: In command ogr2ogr add ogc_fid in output

texts:
I am executing the command

ogr2ogr -nln "demodata" -lco GEOMETRY=AS_XY -t_srs EPSG:4326  -preserve_fid
        -f "CSV" "/tmp/demo.csv" "PG:host=localhost dbname=demo user=root password=root schemas=public tables=demodata" "public.demodata"

I am getting all fields except FID (ogc_fid). That works for my requirements when I create shapefiles, and I don't want to change source elements. But for CSV file I also need the FID.

I tried to add the parameter to the command:

-addfields "ogc_fid"

In both cases I get:

ERROR 1: Field 'ogc_fid' not found in source layer.

Is there any way to force ogr2ogr to return FID no matter if it is in structure or not.


texts:
I am using ArcGIS Desktop 10.5.1.

I have a raster image which I have extracted from a large raster data, and ended up having some scattered small raster groups on the extracted image. This scattered raster are very difficult to deal with when changing the raster to polygon.

How can I delete scattered raster data sets which are of no use?


texts: Mask all objects less than given pixel count

texts:
I'm trying to remove all objects from my analysis that are less than a given size: in this case, less than 36 pixels. I had thought applying image.gt(36) would give me what I need, but when I use that to update the image, I don't see a different result.

How can I access the  object size to mask small objects?

// set the map to our region of interest
var pWestfield = ee.Geometry.Point([-72.87123222351073, 42.29052524043022]);
Map.setCenter(-72.87123222351073, 42.29052524043022, 13);

// get our dataset
var iGlobalSurfaceWater = ee.Image("JRC/GSW1_0/GlobalSurfaceWater");

// create a binary version of occurrence: did it or did it not?
var justOccurence = iGlobalSurfaceWater.select("occurrence");
var binary = justOccurence.and(justOccurence);
Map.addLayer(binary, {palette: 'FF0000'}, "Binary");

// Compute the number of pixels in each patch.
var patchsize = binary.connectedPixelCount(1024, false);
Map.addLayer(patchsize, {palette: ['000000', 'FFFFFF']}, 'patch size');

// Get rid of any patches that are less than 36 pixels.
// Theoretically.
var patchsizeLarge = patchsize.gt(36);
patchsizeLarge = patchsizeLarge.updateMask(patchsizeLarge);
Map.addLayer(patchsize, {palette: ['000000', 'FFFFFF']}, 'patch size large');


texts: Can I view ArcScene 3D models in a free viewer such as ArcReader or another 3rd party viewer?

texts:
I would like to build a complex 3D model in ArcScene, export it and send it to someone else that does not have ArcGIS with 3D Analyst. The idea is to interactively navigate the model during a presentation.

I've never used ArcReader or ArcGIS publisher before so I have no idea if a complex 3D model can be published and exported.

If not what are my alternatives for viewing the completed model outside if ArcGIS?


texts:
I have a question about projections. I have a hillshade raster with UTM 31N projection (I reprojected this from WGS 84). I also have a satellite image that I need to georeference. For georeferencing I choose UTM 31N as projection for the modified satellite image, however when it opens in QGIS, it is not in the same spot as the hillshade raster. When I change the projection of the modified satellite image to WGS 84, it works well and they overlap. What did I do wrong?

UL: 41˚46'40'' / 0˚30'19''
LR: 42˚37'28'' / 1˚01'37''


texts: Extra Column after Using Polyline to Raster

texts: Does anyone know why there is an extra column in the resulting raster file (Raster2) after using the tool Polyline to Raster in ArcMap 10.2? I had imported the cell size and set the extent to be the same as another raster file (Raster1). In the end. Raster 2's cell size and number of rows is the same as Raster 1, but Raster2 has one extra column.


texts: Combining and displaying demographic data?

texts:
I am currently undertaking a demographic research project for my employer. We have the data on income, age, sex, language, transport etc and I have mapped most of it and created some thematics of said variables.

Whilst the maps are perhaps adequate from a raw data representation perspective, I am looking for a technique or thematic style which displays statistical and demographic data in a more intuitive and interpretive fashion.

Here is an example of what i would like to improve:

Is there a specific technique or method that people are using which displays thematic data which is more easily interpreted by a non GIS user?

Similar to this question but perhaps with more technical descriptions on how to achieve some of the better cartographic outputs.

@IanS suggested using rasters to generate a heat map type approach:


texts: Differences and thresholdis of ImageCollection in GEE

texts:
I have the following code (very basic as I am new user), and I want to create an NDVI difference collection by subtracting every NDVI image from my initial collection, with its previous one (something like NDVI2-NDVI1,NDVI3-NDVI2,...).
After that, I want to do some thresholding and assign specific values, but I do not know how to do this.

My attempt:

var imageCollection = ee.ImageCollection("MODIS/006/MCD43A4")
 .filterBounds(geometry)
    .map(function(image){return image.clip(geometry)})
      .filterDate('2005-07-01','2005-08-31');
   var NDVIbands = imageCollection.select("Nadir_Reflectance_Band2","Nadir_Reflectance_Band1","BRDF_Albedo_Band_Mandatory_Quality_Band1","BRDF_Albedo_Band_Mandatory_Quality_Band2");
   var geometry = /* color: #98ff00 */ee.Geometry.Polygon(
        [[[-53.90538056149785, -3.6610951156347706],
          [-53.860061957982225, -3.959810400282152],
          [-53.362930610325975, -3.8529425805852524],
          [-53.42198212399785, -3.5706383493989544]]]);

var masks = function(image){
  var mask1 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band1").eq(0);
  var mask2 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band2").eq(0);
  return image.updateMask(mask1).updateMask(mask2);
};

var filtered = NDVIbands.map(masks);

var myNDVI = function(image) {
  var ndvi = image.normalizedDifference(['Nadir_Reflectance_Band2', 'Nadir_Reflectance_Band1']).rename('NDVI');
  return image.addBands(ndvi);
};

var NDVI = filtered.map(myNDVI);

//From this point, I need your help because I am not familiar with this.

//I do not know how to set the slices in order to subtract every NDVI from its previous
//(like NDVI2-NDVI1,NDVI3-NDVI2,...).Based on what I want, my NDVIdif collection must have
//one element less than the initial NDVI collection.

var NDVIdif = NDVI.slice(,).subtract(NDVI.slice(,))

//And here, I want to create classes based on value-limits but again, I do not know how to do this.

var reclass = function(img) {
  var mclass = img.remap([-0.250 - -0.103,-0.105 - 0.201],[  1,   2],   null, 'reclass dNDVI');
  return mclass;
};


texts: Styles not applied when using QgsMapRenderer

texts:
Using pyqgis, i am trying to export some maps as images using QgsMapRenderer. Here is the layer/style setup, both load correctly:

lyr = QgsVectorLayer(uri.uri(), 'Report', 'postgres')
styleload = lyr.loadNamedStyle(stylepath)
if not styleload[1]:
    print 'Error loading report selected style'

if not lyr.isValid():
    raise IOError, "Failed to open the reports layer"
else:
    QgsMapLayerRegistry.instance().addMapLayer(lyr)

Here is how i export, using QgsMapRenderer...

renderer = QgsMapRenderer()
renderer.setDestinationCrs(QgsCoordinateReferenceSystem(crs))
renderer.setLabelingEngine(QgsPalLabeling())
renderer.setLayerSet(QgsMapLayerRegistry.instance().mapLayers().keys())

p = QPainter()
extents = lyr.extent()
extents.scale(1.1)
renderer.setExtent(extents)

img = QImage(QSize(2050, 1150), QImage.Format_ARGB32_Premultiplied)
img.fill(QColor(255, 255, 255).rgb())

p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

renderer.setOutputSize(img.size(), img.logicalDpiX())
renderer.render(p)

img.save(mappath, 'jpg')
p.end()

The problem is, i get the image out, but styles are not being applied at all, even though they load correctly.


texts: QGIS swaps(?) attribute tables of PostGIS layers

texts:
I have a QGIS project (2.18.14 Win10), last saved yesterday, with several PostGIS layers. As I opened the project today the PostGIS layers did not display on the map (but appear in the Layers Panel). The layers' properties (source, etc.) are correct, but their attribute tables are those of other PostGIS layers. It's as if QGIS somehow mixed up the attribute tables.

Other layers (shapefiles) in the same project display correctly.

I have not touched the database. An ArcMap document saved at the same time and accessing the same PostGIS layers shows no issues.

Any ideas on what might be going on and how to fix this?


texts: Investigating PostGIS table not showing in QGIS

texts:
I created a table in Postgres:

create table eric_zs as
select z.id, z.name, st_multi(st_union(a.geom)) as geom, count(*)
from geo_as a
     inner join geo_zs z on a."ZS_id" = z.id
group by z.id, z.name;

That table is showing fine in QGIS.

Now, if I replace the st_multi(st_union()) with for example ST_Collect(ST_ExteriorRing(a.geom)), the table just won't show up in QGIS.

My question is about how I should go about investigating this kind of issue. Is there a log, a specific condition to check to get a clue about why it is not showing up ?


texts:
Under Settings>Options>Layout>Layout Paths we can define Path(s) to search for extra print templates like shown here:

I intend to deploy an organization-wide standard path for print templates using PyQGIS in a startup-script but struggle in finding a method that does what I desire.

I cannot find something related in the QgsSettings(), but I found a method

QgsApplication.instance().layoutTemplatePaths()

returning a list of the currently defined print template paths, but no equivalent method setLayoutTemplatePaths().

How may I set these print template paths using PyQGIS?


texts: Create/deploy REST API using ArcGIS/Python

texts:
I've written a bit of Python code that utilizes ArcGIS to perform some analysis on a specific location. Running this locally on my computer works great. However, I am struggling and am unfamiliar with how to deploy the code as a REST API that I can access through a website.

What is the best or easiest tool to get this done? I've read about some more complex features of ArcGIS such Enterprise version which costs tens of thousands of dollars a year, but that hardly seems reasonable for a tiny one man passion project.

High level I would want the process to be:
1) User inputs address into website using REST API
2) GIS & Python code runs (utilizing ArcGIS scripting tools and some data)
3) Returns to user the results of the analysis through the website.


texts: Leaflet's worldCopyJump breaks map, zooming fixes it up (updated)

texts:
I had my map working with worldCopyJump set to true, but then it stopped working. I'm not sure exactly what broke it (the original form of this question said there was a connection to a plugin, but that doesn't seem to be the case), but it might be that I updated to Leaflet 0.62.

In any case, this is the problem: when the map loads, dragging the map would throw a bunch of errors claiming to have invalid LatLong objects (with a NaN instead of the longitude), and the map would freeze.

Oddly enough, if before dragging the map I zoomed in or out, and then dragged, everything would work fine. Whatever the map is doing during onzoom seems to be somehow fixing whatever problem it had.

Also peculiar is the LatLongs that are reported to be invalid: they are nowhere in the code. I have no idea where those numbers might be coming from.

I traced the error back to the worldCopyJump option. If this is set to false then the map works from the get go. But I need that option turned on. Can anybody think of what might be going on?

You can find a (non-)functional version of the map online here.

In case that does not last, here's the code:

  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.6.2/leaflet.css" />


texts: Calling ST_Transform in geometry that is already in target srid?

texts:
I'm writing a function that has a geometry as a parameter that can be in different projections:

CREATE OR REPLACE FUNCTION foo(in geometry) RETURNS boolean AS $$
DECLARE
    transformed_geom geometry := in;
BEGIN
    IF ST_SRID(transformed_geom) != 32737 THEN
        transformed_geom := ST_Transform(transformed_geom, 32737);
    END IF;
    [ ... ]
END;
$$ language plpgsql;

The function needs the geometry in a specific projection (32737) but I don't want to obligate the caller to make the transformation, I prefer to handle it inside the function.

The question is, instead of check the actual srid, it's safe to always call ST_Transform? Will do anything if the geometry is already in the target projection?

I check with a query like this that it seems to be safe, but it will be awesome if someone can confirm this as I did not see it in the docs.

SELECT ST_AsEWKT(geom), ST_AsEWKT(ST_Transform(geom, 32737)) FROM table_in_32737;


texts: How to limit calculations of r.grow.distance to contiguous cells in raster mask?

texts:
I have a layer with rasterized vector polygons (white cells) in a GRASS database. The region is set to the bounding box of the polygons. I also set a raster mask to mask all cells between these polygons (grey cells). In addition, I have a second layer with rasterized linestrings which a located inside the polygons (coloured cells).
Now I want to use r.grow.distance to calculate the distance of every raster cell to the next raster feature BUT only inside the same contiguous "zone". For example, raster cells of polygon 1 should only be related to the rose raster feature. Also, polygon to would be assigned to the purple feature, although there is no raster feature on it. Currently, the south-eastern cells of polygon 1 would be associated with the purple raster feature, because it's closer than the rose on.
How can I force r.grow.distance to limit the calculation to each contiguous zone and not to "jump" across masked cells?
I would like to avoid running this command for every polygon separately as a possible work-around.


texts:
I created a LAS Dataset out of my 1000 LAS-Files (regular grids) with the purpose of converting the LAS Dataset to a Raster later on.

But the problem is: my LAS-Files are overlapped on the borders and I want the LAS Dataset to maintain all the values at the borders so I can average them later on using "LAS Dataset to Raster". Unfortunately my created LAS Dataset only chose one value to keep.

Anyone has an idea of how to keep overlapped values on the LAS Dataset?


texts:
I have a Postgres 8.4 PostGIS 2.0 enabled database and I've loaded it up with a bunch of shapefiles.

I'm trying to get the database to do on-the-fly unions of geometries.

Here is the query I'm working with:

SELECT ST_AsGeoJSON(ST_Union(the_geom)) FROM (SELECT geom FROM forecast_zones WHERE
state = 'ND' AND zone IN ('002','011','021','035','046','050')) as the_geom;

The Subquery to pull the geometries of the individual zones works fine:

SELECT geom FROM forecast_zones WHERE state = 'ND' AND zone IN
('002','011','021','035','046','050');

But whenever I try to pass that to ST_Union I get undefined function errors.

EDIT:

Here is the query I ended up using that works:

SELECT ST_AsGeoJSON(ST_Union(forecast_zones.geom)) AS geom FROM forecast_zones
WHERE forecast_zones.state = 'ND' AND forecast_zones.zone IN
('002','011','021','035','046','050')
GROUP BY forecast_zones.state


texts:
Are there any free global address datasets I can use to test my geocoding application?  I'm looking for a dataset that contains residential addresses from countries around the world.  I don't need a name, just a street address, region, country, and postal code (if the country uses postal codes).  If it has a latitude and longitude for each address, that would be a plus.

I saw this question but I'm wondering if anything has changed since it was answered.

Thanks.


texts: What are the most efficient Geographic Weighted Regression packages?

texts:
What are the most memory efficient open source packages for calculating a geographically weighted regression (GWR)?

I am in a situation where I need to do a geographically weighted regression on a set of points where training data consists of about 40,000 observations and each observation has about 20,000 variables.

I have implemented a version of GWR myself using a combination of Python Numpy/SciPy and PostGIS.  I solve the regression using a matrix algebra approach, but this fails due to memory issues when I have dense, feature rich systems with many observations.

One way to get around the memory issue is use an iterative approach for finding a line of best fit, such as an incremental gradient descent.  I'm thinking it should work something like (http://www.eecs.wsu.edu/~cook/dm/lectures/l5/node14.html).  Incremental Gradient Descent is described pretty well here in pages 4-7 (http://cs229.stanford.edu/notes/cs229-notes1.pdf).

Obviously I could implement this myself, but I was hoping maybe someone else had already coded something similar.


texts: Adding unique ID for groups of items using field calculator in QGIS

texts:
I have a set of string values in a shapefile's attribute column and want to assign a number to them in a way that each unique string gets an integer number, like this:

How do I do this using the field calculator? I have 3,000 rows of data and 97 different names to do this for.


texts: Compare distributions of point patterns for significant differences in clustering

texts:
My problem is I have tens of thousands of points (GPS) representing animal locations over a period of several months. I want to test if a portion of those points (only about a few hundred) which occurred during stressful events have a markedly different distribution than the overall distribution of the remainder points.

What method/statistical route should I use? I've heard of Ripley's K Function, but I don't know enough about it to know if it's applicable, especially with the relatively small amount of points representing the disturbed case. The trouble is also that many of these methods sound like they indicate if a distribution is clustered or random but not whether or not distributions seem to match.
For this test I would have time involved but I will most likely assume that their behavior is not seasonality-dependent (although I know inevitably it is)
If it's possible I would also like to compare whether the directions the animals tended to move varied between the two circumstances. So normally the animals follow a particular pathway and locations off of the pathway tended to result in motion towards the pathway, I want to know if I can test to see if the disturbed points tended to have their directionality (as determined by their next subsequent point) no longer fit that pattern, i.e. avoiding the pathway.


texts:
Does ArcGIS online have the capability of providing real time (live) tracking of assets?  If I have a staff member with a cell phone can ArcGIS online show me where that person is at any given moment?

I have ArcGIS desktop.  I do not have ArcGIS Server.


texts: Free download NDVI images or raw satellite images for Sri Lanka?

texts: Could anyone provide a site to download NDVI images or raw satellite images of Sri Lanka (Asia) for free for me to calculate the NDVI.


texts: Differences in MaxEnt between 'Biomod2' and 'Dismo' package

texts:
I am trying to run an ensemble species distribution, however the the results I am getting for MAXENT.Phillips model are not reasonable projections compared to the Random Forest and Boosted Regression models that I am using.

The projection for MAXENT deviates heavily from the other projections. I attached a picture demonstrating first the Random Forest projection, then the Boosted Regression, and on the right is the MAXENT projection, which appears to not be strict at all with it's predictions.

[

Using the same environmental variables, presence and background points, I ran a projection in the 'Dismo' package and obtained results more reasonable and in agreement to the Random Forest and Boosted models in BIOMOD2.

This leads me to believe that I have not adjusted a default setting in BIOMOD, but I am at a loss for what I need to do.

I included my code below. Why am I getting such different projections in MAXENT?

For the environmental variables, I used the getdata function.

ext.1<-extent(95,115,25,45)
bio.15 <- getData('worldclim',var='bio',res=2.5)
bio.15<-crop(bio.15,ext.1)
bio.15<-stack(bio.15)

And for the projection environmental variables:

bio.70<-getData('CMIP5',var='bio',res=2.5,rcp=85,model='AC',year=70)
bio.70<-crop(bio.70,ext.1)
bio.70<-stack(bio.70)
names(bio.70)<-names(bio.15)

Beyond that, I kept everything relatively similar to the BIOMOD2 vignette. The link to the csv is here csv_link:

file<-paste([csv_link])
data.short<-read.table(file,header=TRUE,sep=",")
resp.name<-'GiantPanda'
resp.test<-as.numeric(panda.data.short[,1])
resp.xy.test<-data.frame(panda.data.short[,c("POINT_X","POINT_Y")])
ensemble.test<-BIOMOD_FormatingData(resp.var=resp.test,
                               expl.var=bio.15,
                               resp.xy=resp.xy.test,
                               resp.name=resp.name)

The .jar file is in the working directory, and there are no spaces in the working directory

ensembleoption<-BIOMOD_ModelingOptions()

ensemble.out.test<-BIOMOD_Modeling(
  ensemble.test,
  models=c('RF','GBM','MAXENT.Phillips'),
  models.options=ensembleoption,
  NbRunEval=1,
  DataSplit=70,
  Prevalence=NULL,
  VarImport=1,
  models.eval.meth=c('TSS','ROC'),
  SaveObj=TRUE,
  rescal.all.models=FALSE,
  do.full.models=FALSE,
  modeling.id="3_28")

Then the projection is:

proj.future<-BIOMOD_Projection(
  modeling.output=ensemble.out.test,
  new.env=bio.70,
  proj.name='projection_future_3_28',
  selected.models='all',
  binary.meth=NULL,
  compress='xz',
  build.clamping.mask=T,
  output.format='.grd')


texts:
I've had an issue with the Processing plugin over that last few days. QGIS
has been really slow to start (about 15 minutes in some cases) and I then
get a message that says "Plugin Processing: the plugin will be disabled
because it crashed QGIS during last startup.

Please report and issue and
re-enable the plugin when the problem has been resolved"

Can anyone tell me how I report the issue? What logs I need, where do I
find them and who do I send them to?


texts: How to install GDAL version 3.0 on Windows in order to generate GeoPDF with Qgis 3.10?

texts:
I am very enthusiastic to create some GeoPDFs using the new version of Qgis 3.10 on Windows. Unfortunately, using the OSGeo Installer, the newest gdal package version proposed is 2.4.1, which is annoying considering that the minimal required version to create GeoPDFs is the 3.0.

As described in the documentation here https://gdal.org/download.html#current-releases, I tried to install the latest gdal version through Anaconda client. Here again the latest package available is 2.3.3.

In googling on the web I also found that site => http://www.gisinternals.com/development.php which gives some links to download gdal 3.0. Despite the fact that the installation is possible that way (I guess), I do not get it. Besides, how could we be sure that the gdal installation that way will be "clean" and correctly working in QGIS 3.10?

I really think that new feature of creation GeoPDFs is a real improvement. I guess I will not be an isolate case, and this is really too bad that it does not work "out-of-the-box". Maybe I am also too eager to use that function, and we have to wait a few dasy/weeks to have the new gdal 3.0 automatically suggested in the OSGeo Installer...

Could someone explain the way to install this gdal 3.0 on Windows in order to create GeoPDFs with Qgis 3.10?


texts:
I am using PostgreSQL 9.3.3 & PostGIS 2.1.1 and I have a table Entries which has a column called "latlng" of geographic type (this column stores entry's latitude and longitude coordinates).

I am trying to find all entries that belong to a given polygon - like this:

entries = Entry.find_by_sql("
    SELECT *
    FROM Entries
    WHERE ST_Within(
        Entries.latlng::geometry,
        ST_GeomFromText(
            'POLYGON((
                #{swLat} #{swLng},
                #{swLat} #{neLng},
                #{neLat} #{neLng},
                #{neLat} #{swLng},
                #{swLat} #{swLng}))',
            4326)
    );")

Which translates into the following SQL statement:

SELECT *
FROM Entries
WHERE ST_Within(
    Entries.latlng::geometry,
    ST_GeomFromText(
        'POLYGON((
            45.61490708797995 -121.32521530382634,
            45.61490708797995 -121.32338469617367,
            45.615492910571655 -121.32338469617367,
            45.615492910571655 -121.32521530382634,
            45.61490708797995 -121.32521530382634))',
        4326)

Even though I am sure that my table has entries that have coordinates that belong to such polygon, I get zero results back. Any ideas on what could be causing this?


texts: I collected data along two perpendicular transects (in a "+" shape) and would like to create a prediction surface map with kriging. What form of kriging (method, parameters, etc.) is appropriate to better model this particular spatial distribution of the data?


texts: OpenLayers clear path on marker movement

texts:
I'm using OpenLayers 6.
I would like to produce behaviour like on Google Maps:
When the marker has passed the point, all previous points should be hidden and the only route that is left to pass should be visible.
For example:
https://openlayers.org/en/latest/examples/feature-move-animation.html
So when I start marker animation, I would like the previous yellow route that marker passed to be hidden, and only route ahead of him should be visible
Is that even possible?


texts: Determine the UTM zone/hemisphere in a particular datum

texts:
I'm using GDAL to convert coordinates from one system (geographic or projected) to another.

Assume I want to convert geographic coordinates (latitude and longitude) in the WGS84 datum, to UTM coordinates in the Amersfoort datum (EPSG:4289).

In order for GDAL to determine the UTM coordinates (easting and northing), it needs:

the input spatial reference system: geographic, based on the WGS84 datum
the input coordinates: latitude and longitude
the target spatial reference system: projected - UTM -, based on the local Amersfoort datum

To instantiate the target SRS, besides the datum, I need to specify the zone and the hemisphere.

The question is: How to I determine the UTM zone and hemisphere where the input geographic coordinates fall in? I could use the Wikipedia article for that. As far as I understood, the mechanism described there works for determining the zone and hemisphere for geographic coordinates in the WGS84 datum. What is not clear to me is whether the resulting zone and hemisphere make sense only for the WGS84 datum, or if they can be used in any, local or global, datum.

If the zone and hemisphere determined from geographic coordinates in the WGS84 datum only make sense for determining the UTM coordinates in that same WGS84 datum, how can I determine the UTM zone and hemisphere in another datum - WGS72, NAD27, NAD83 etc.?

Also, does it even make sense to speak of UTM coordinates in a local datum?

I hope I haven't messed up the concepts too much. Anyway, any pointers will be greatly appreciated.


texts: Recurring system crash while using GDAL Georeferencer in QGIS 2.18

texts:
I am currently working on aerial images I have to reference in QGIS. I usually use the GDAL tool for that, but during the last days, QGIS crashes all few minutes whenever I use it. I ran into this problem before, but usually a restart of the computer fixed it for some reason. This time, that did not help.

Does anyone know this error and how to avoid it? This problem really gets onto my nerves.

I am using QGIS 2.18.17 on a 64 bit Windows 7 PC.


texts:
I have a code in which I'd like a da.UpdateCursor, I have probles with the "where clause", here I am about to use a variable that I get from arcpy.Describe. I have no idea of this even possible.
My code is:

    newValue =  arcpy.GetParameterAsText(0)
    inlyrdsc = arcpy.Describe("layertoedit")
    inFID = inlyrdsc.FIDset
    num = inFID.split("; ")

    inlayer = arcpy.mapping.Layer("layertoedit")
    query = "[OBJECTID] = %s", (num)
    fieldName = "fieldtoupdate"
    with arcpy.da.UpdateCursor(inlayer, (fieldName,), query) as cursor:
        for row in cursor:
            row[0] = newValue
            cursor.updateRow(row)

I get the following error message:

TypeError: 'where_clause' is not a string

Anyways my goal is to write a script which fills the attribute table of the selected element (the selection is described by the OBJECTID).
This might not be the best solution so any ideas are welcomed.


texts: For Each loop over an Array in MapBasic?

texts:
.NET C# developer here. Trying to figure out MapBasic with the IDE for the first time. It's shall we say...a different experience to C# in Visual studio.

I have a ButtonPad with a button that calls a stored proc in the SQL server db which performs an insert on one selected polygon in my BidBlock layer:

' User Hit Bid Block
print cmd_INFO_MENUITEM
if getSelectionInfo(SelCount, SelTableName) then
    if SelCount = 1 then
        'Have one Contract selected
        if selTableName = vwBidBlocks then
            i_blockId = selection.block_id
            if i_blockId = 0 then
                run menu command 304
                note "Please re-select your Bid Block."
        '       set handler WinChangedHandler on
                exit sub
            end if
            'call notyet()
            call BidBlock(i_blockID)
        else
            dialog error handling blah blah

I'm now trying to implement this function to perform on a selection of > 1 polygons in the layer. I want to therefore make an dynamic array of items using reDim(selection)?

& itterate over the ID's of each selected item in the array calling my BidBlock(i_blockId) function. In C# I would use a for each loop here which implements the IEnumberable interface. Is there any way of doing this in MapBasic?


texts: Export MapInfo map to Google Earth (.kml file)

texts:
I have a map in MapInfo.

I can convert it to Google Earth file (.kml) manually but I would like to do it automatically by commands in a function.

Can anyone help me?


texts: Offset in GeoTIFF coordinates Google Earth Engine extracts

texts:
I download a certain geographic location (lon,lat) with buffer radius around for MODIS products, such as MCD43A2/MCD43A4 from Google Earth Engine using ee, such as:

 lon = 10.452175
 lat = 51.079206
 product = 'MODIS/006/MCD43A4'
 sensor = product.split('/')[-1]
 radius = 2000
 res = 500
 start,end = (datetime(2004, 1, 1, 0, 0), datetime(2004, 12, 31, 23, 59))

geometry = ee.Geometry.Point(lon, lat).buffer(radius)
col      = ee.ImageCollection(product)
features = col.filterDate(start, end).getInfo()['features']

http = urllib3.PoolManager(num_pools=10,
                            cert_reqs='CERT_REQUIRED',
                            ca_certs=certifi.where())

for feature in features:
    url = ee.Image(feature['id']).clip(geometry).getDownloadURL()
    #request.urlretrieve(url,OutDir.format(site_id,sensor)+feature['id'].split('/')[-1]+'.zip')
    r = http.request('GET', url)
    with open(OutDir.format(site_id,sensor)+feature['id'].split('/')[-1]+'.zip', "wb") as outfile:
        outfile.write(r.data)

than I get, using gdalinfo 2004_01_06.Nadir_Reflectance_Band1.tif :

Driver: GTiff/GeoTIFF
Files: 2004_01_06.Nadir_Reflectance_Band1.tif
Size is 10, 10
Coordinate System is:
PROJCS["MODIS Sinusoidal",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (727864.277675088495016,5682067.155503926798701)
Pixel Size = (463.312716528000010,-463.312716527000021)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  727864.278, 5682067.156) ( 10d25'45.53"E, 51d16'10.18"N)
Lower Left  (  727864.278, 5677434.028) ( 10d25'11.63"E, 51d13'40.26"N)
Upper Right (  732497.405, 5682067.156) ( 10d29'44.52"E, 51d16'10.18"N)
Lower Right (  732497.405, 5677434.028) ( 10d29'10.41"E, 51d13'40.26"N)
Center      (  730180.841, 5679750.592) ( 10d27'28.01"E, 51d14'55.22"N)
Band 1 Block=10x10 Type=Int16, ColorInterp=Gray

This does not refer to the input I made:
'lon = 10.452175, lat = 51.079206'
at least for latitude.

How can I get around this, retrieving the original (projection) data subset, not re-projected, with the center (lon,lat) in the center of the subset? Or which would be the correct transformation to e done?


texts:
I work with arcview 10.3. Every time i export maps with ecw raster format , python (version 2.7.8) get crashed and i see:

I use this code:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\yaron\shonot\software\gis\maps\gis"
counter = 0
for mxd in arcpy.ListFiles("*.mxd"):
    print (mxd) # print list of mxd's in the folder
    mapdoc = arcpy.mapping.MapDocument(r"C:\yaron\shonot\software\gis\maps\gis\\" + mxd)
    arcpy.mapping.ExportToJPEG(mapdoc, r"C:\yaron\shonot\software\gis\maps\gis\output\\"+os.path.basename(mapdoc.filePath).replace ('.mxd','') +".jpg",resolution=196)
    print ('ExportToJPEG')
    mapdoc.save()
    counter = counter + 1
del mxd
print (counter)


texts: Getting co-ordinates of current map canvas using PyQGIS?

texts:
I'm trying get the extents of the current map canvas. CRS is 22700 BNG.
My code:
canvas = iface.mapCanvas()
print canvas

vLayer = iface.activeLayer()
extents = vLayer.extent()

print extents

Which returns

qgis._gui.QgsMapCanvas object at 0x1198c6a68
qgis._gui.QgsMapCanvas object at 0x1198c6a68

I'm want 4 numbers: xMin, xMax, yMin, yMax


texts: Cause and effect of a large number of 'not in feature positions' in trimble .cor files

texts:
I have been collecting GPS data in the field using a Trimble G series nomad. The collection site is a forested area with moderate tree cover. The features will be exported into a shapefile for further analysis in ArcGIS. However, Following differential correction some of my features are still above the 1m horizontal accuracy I require so I am am manually removing the feature related positions with very high vertical accuracy in an attempt to reduce the effect of multipath errors.

Upon opening the files in PFO I am noticing a large number of not in feature positions (as seen in the attached picture) which are complicating the process as they overlap with the feature related positions upon feature deletion, so I am wondering:

If the these positions are indicative of a methodological misstep in my collection?
what effect, if any these may have on the horizontal accuracy of the features and as such if deleting them is likely to negatively impact the accuracy of the features?


texts: Creating raster by randomly choosing cell value from multiple overlapping rasters?

texts:
I'm using ArcGIS Desktop 10 with its Spatial Analyst extension.

How do I combine multiple rasters into one, always choosing randomly from the values of overlapping cells?

I have an image that may explain this better:


texts: ArcGIS Pro calculate field using arcade not working?

texts:
Im trying to calculate a field called "bark t/ha" using the following expression:

if ($feature.Barkfuelhaz =="Low"){
  return"0"}
else if($feature.Barkfuelhaz =="Moderate"){
  return"1"}
else if($feature.Barkfuelhaz == "High"){
  return"2"}
else if($feature.Barkfuelhaz == "Very High"){
  return"5"}
else if($feature.Barkfuelhaz == "Extreme"){
  return"7"}

the expression validates however when I run it, all values in the Bark/ta field change to "0"

I've even tried running a simpler one "if ($feature.Barkfuelhaz == "Moderate"){
  Return "2"}"
but all values still stay as 0


texts: Calling a function after map display rather than map load

texts:
I have the following method which I am calling on load of the map. However when I try to turn on a particular layer within this code I am unsuccessful as although the map would have been loaded it would not be displayed. Is there some way to call the method after the map is visible on the display? I also tried including the method in document.ready() however it was also unsuccessful.

      app.map.on("load", function () {

                var coords, zoomLevel;
                var urlObject = esri.urlToObject(document.location.href);

                if (urlObject.query && urlObject.query.coords && urlObject.query.zoomLevel && urlObject.query.LotNumber && urlObject.query.LayerName) {
                    var coords = urlObject.query.coords.split(',');

                    var lon = parseFloat(coords[0]);
                    var lat = parseFloat(coords[1]);
                    var lotNumber = urlObject.query.LotNumber;
                    var layerName = urlObject.query.LayerName;

                    var zoomLevel = parseInt(urlObject.query.zoomLevel);
                    var point = new esri.geometry.Point(lon, lat, new SpatialReference({ wkid: 23033 }));

                    for (var i = 0; i < app.map.graphicsLayerIds.length; i++) {
                        var tempLayer = app.map.getLayer(app.map.graphicsLayerIds[i]);

                        var lname = layerName.split(".");
                        layerName = lname[lname.length - 1];

                        if (tempLayer.url.indexOf(layerName) != -1) {
                            tempLayer.show();
                        }
                    }

                    app.map.centerAndZoom(point, zoomLevel);
                }
         });


texts: How to convert .shp to lat and long coordinates en excel?

texts: I am trying to convert the below SHP file into lat and long coordinates in excel. I have no idea how to code or even find where the coding section is.


texts:
I've got plans of my house and a warehouse on paper. I would like to scan them in JPEG, import them on my PC and measure them.
I'm looking for a website where I can upload those tiff, jpeg or pdf files, define a scale (measuring something of known measure to give a correct scale to the entire image), measure the rooms by drawing polygons over them, and print out the images with my measured polygons or send them to my email with the total area reported.
I can't find anything like this on the web. Is this kind of utility uncommon? Can any body suggest a simple solution for measuring my planimetry without installing any software?


texts: How to detect overlapping features from 2 different layers with Mapbox GL JS

texts:
I have 2 layers in total,

one is a point layer where each point represent the location and
attributes of a company,
the other is a draggable polygon layer which is regarded as a
filter for spatial analysis for point layer    (with reference from
Mapbox example and here)

The orange points below are companies and the blue circle is the draggable filter area.
What I am going to do is to allow users to drag the blue circle around, so they can understand what companies (and their attributes), are included under the area covered by the circle.
Is Mapbox GL JS able to get the points that are overlapping with the blue circle?


texts: XYTools truncating decimals on export to excel

texts:
Is there some sort of setting to prevent XYTools from truncate decimals when exporting a shapefile to excel?

I live in Sweden hand we normally use comma as decimal point. Even if I change the setting both on system level as in Excel - XYTools still truncate decimals.