texts: Automate Selection/Copy Features from Selection

texts:
I am new to Python and ArcPy, but I would like to automate this process and/or make my code more efficient here. I am basically selecting a layer by attribute, then copying features from the (selected) layer, and then clearing the selection and starting over again with another layer. I need to do this for many layers in a row in my map document. I have attached my Python workflow for two layers below as an example.

>>> arcpy.SelectLayerByAttribute_management("Outlet_A", "NEW_SELECTION", "Struc_Rec = 'Major Repair'")
<Result 'Quad A\\OM_Rec\\Outlet_A'>
>>> arcpy.CopyFeatures_management("Outlet_A", "MajorRepairOutlet_A")
<Result 'X:/data.gdb\\MajorRepairOutlet_A'>
>>> arcpy.SelectLayerByAttribute_management("Outlet_A", "CLEAR_SELECTION")
<Result 'Quad A\\OM_Rec\\Outlet_A'>
>>> arcpy.SelectLayerByAttribute_management("Inlet_A", "NEW_SELECTION", "Struc_Rec = 'Major Repair'")
<Result 'Quad A\\OM_Rec\\Inlet_A'>
>>> arcpy.CopyFeatures_management("Inlet_A", "MajorRepairInlet_A")
<Result 'X:/data.gdb\\MajorRepairInlet_A'>
>>> arcpy.SelectLayerByAttribute_management("Inlet_A", "CLEAR_SELECTION")
<Result 'Quad A\\OM_Rec\\Inlet_A'>


texts:
First time poster so I hope I did this right.
I'm relatively new to Python scripting and am trying to understand the error I'm getting.

I have the following simple code and what I don't understand is why PAD_HI_250extent generates the Method fields does not exist error?  They are both polygon feature classes stored in the same file gdb.  Is there a difference in the two feature classes that I'm not aware of?  Can the content of the fields, as opposed to the field name themselves, affect this (I do know that there are special characters, e.g.,  K'au Forest Reserve)

Running ArcMap 10.1 on Win7, Python 2.7.2

=================================================================================
SCRIPT
======
import arcpy, time

## set environment settings
arcpy.env.overwriteOutput=True
arcpy.env.workspace = r"F:\USER\FLATHER\TE_AtRisk_grid_model\FeatureClasses.gdb"

files = ["BOC_HI_BG_250extent", "PAD_HI_250extent"]

for file in files:
    print 'File:', file
    desc = arcpy.Describe(file)
    fieldnames = [f.name for f in desc.fields if f.type not in ["Geometry", "Raster", "Blob"]]
    print fieldnames
    print ' '

==================================================================================
OUTPUT
======
File: BOC_HI_BG_250extent
[u'OBJECTID', u'FID_HI_fishnet_250sqmi', u'FID_BOC_HI_BlockGroups', u'STATEFP', u'COUNTYFP', u'TRACTCE', u'BLKGRPCE', u'GEOID', u'NAMELSAD', u'MTFCC', u'FUNCSTAT', u'ALAND', u'AWATER', u'INTPTLAT', u'INTPTLON', u'Shape_Length', u'Shape_Area', u'tot_area', u'fishnetID']

File: PAD_HI_250extent

Traceback (most recent call last):
  File "F:/USER/FLATHER/TE_AtRisk_grid_model/scripts/DDD_descfields.py", line 15, in <module>
    fieldnames = [f.name for f in desc.fields if f.type not in ["Geometry", "Raster", "Blob"]]
AttributeError: DescribeData: Method fields does not exist
>>>


texts: How to name a lot of point of a KML file?

texts: I have a KML file with 5400 points with no name. I would like to name each point. My idea was to convert this KML file to CSV to rename quickly them and to convert back this CSV to KML, but I don't succeed doing it.


texts:
I'm trying to Georeference an image. I've followed the guide from QGIS, but for some reason, it's stretching, and offsetting the image from the points.

For example, my input of the southwest corner is:

49.729628
98.471569

And after I've run the georeferencer the southwest corner is:

49.7312
98.9994

Both the image, and the georeferencer are using WGS 84 EPSG:4326 as the CRS.

You can see the input, the settings and output below:


texts: How do I create midpoints for a set of polylines in QGIS?

texts:
I am trying to create midpoints for a set of polylines in QGIS. My python knowledge is poor and was wondering if there was a plugin that did if for me.

I have tried using the Mean Coordinates function suggested in Creating point to center of polyline? but this returns a centroid which is sometimes off the line (especially when lines are not straight). I am essentially looking for a point that is on the line half way along.


texts:
This is a special case, I have many popups attached to objects and they work fine and I am having a problem in only one case.

All layers and popups are handled by select control's onFeatureSelect, but this one is handled by layers featureSelected event.

On opening popup it throws:

  Uncaught TypeError: Cannot read property 'lat' of null.

Expanding error shows error on function OpenLayers.Bounds.OpenLayers.Class.determineQuadrant but the popup opens.

On Clicking the close button, the click event is also fired and onPopupClose function is also called ( I have checked it with monitor and monitorEvents in console) but the popup does not close.  The popup in other layers is still function normally.

Here is the code:

select_Control = new OpenLayers.Control.SelectFeature([business, schools, healthfacilities, open_space], {
    onSelect: onFeatureSelect,
    onUnselect: onFeatureUnselect,
    click: true
});

open_space.events.on({
    "featureselected": function(e) {
        var PopupPos1 = e.feature.geometry.getBounds().getCenterLonLat();
        popup_content = '

sifal ground

';
        var popup1 = new OpenLayers.Popup.FramedCloud("chicken",
            PopupPos1,
            null,
            popup_content,
            null, true, onPopupClose);

        e.feature.popup = popup1;
        map.addPopup(popup1);
        popup1.setBackgroundColor('#FFBBBB');
        popup1.updateSize();
        popup1.draw();
    }
})

function onPopupClose(e) {
    select_Control.unselectAll();
}

function onFeatureSelect(feature) {
    //debugger;
    var PopupPos1 = new OpenLayers.LonLat(feature.geometry.x, feature.geometry.y)
    popup_content = feature.attributes.name;
    var popup1 = new OpenLayers.Popup.FramedCloud("chicken",
        PopupPos1,
        null,
        popup_content,
        null, true, onPopupClose);
    feature.popup = popup1;
    map.addPopup(popup1);
    popup1.setBackgroundColor('#FFBBBB');
    popup1.draw();
}

function onFeatureUnselect(feature) {
    map.removePopup(feature.popup);
    select_Control.unselectAll();
}


texts: Publishing feature layer with enabled editor tracking overwrites edited creator attribute

texts:
trying to publish a populated feature layer via ArcMap 10.2.2 (share as > service) to ArcGis Online with Editor tracking enabled, without having all "Creator" and "Editor" Attributes being overwritten afterwards.

As asking questions in such places is a new thing for me I attached 2 pics to demonstrate the problem:


texts:
In the layer properties of each layer there is an option "Show MapTips using the display expression" that is by default disabled and I would like to find out if there is a way to enable this option by default for all layers or by using arcpy.
I have tried to run arcpy.UpdateLayer but that would update either only the symbology properties or every layer property.


texts: How can you refresh the user attributes in a FME Desktop workspace?

texts:
How can you refresh the user attributes in a workspace?

I have a FME workspace that is reading from a MapInfo table, and writing to SQL Server spatial table

The SQL table has had a new field added to it.

I want the new field to appear in the list of user attributes, but I can't see any "refresh" options. I can't believe it's a case of deleting the existing SQL writer and adding a new one


texts:
I am running the shortest_path function but it will not return rows of lines to follow from source id to target id. Is there a specific way to digitise my dataset? Does pgrouting understand intersections between lines? I digitised my lines one by one connecting snapping them together.

To be a bit more clear, in my db i have a gid for every line, plus a source and target id for that line (plus the other columns). Shouldn't the shortest path function take me from a random source id to a random target id?

I hope my problem is clear but let me know if further explanation is needed :)

Thanks in advance.
A


texts: How to change the settings for where ArcGIS splits world data

texts:
I have a world map made up of 6 shp files.  All the files and the data frame are set to WGS84.  The way Arc 10 and 10.1 decides to display the files, it splits Russia into 2.  The majority of Russia can be seen to the East, however, one small bit is displayed at the western extent of the map.
This is a single polygon.

Is there a way to change the way ArcGIS treats this data so the the whole of Russia is displayed to the east? (I will eventually need to print the data).

EDIT:  I'm starting to think that in order to do this I might need to create a new coordinate system (a copy of wgs84) and offset the centre of it so that it sits to the east.  Thoughts?


texts: Geoserver 2.10 disk quota subsystem disabled

texts:
Since a couple of days, I get this error in my Windows Geoserver 2.10 installation:

Could not open JDBC Connection for transaction; nested exception is org.apache.commons.dbcp.SQLNestedException: Cannot create PoolableConnectionFactory (Unique index or primary key violation: PRIMARY_KEY_B ON PUBLIC.TILEPAGE(KEY); SQL statement: CREATE PRIMARY KEY ON PUBLIC.TILEPAGE(KEY) [23001-119])

Has anyone an idea how to solve this?


texts: TIFF image from R displaying only white in ArcMap?

texts:
I cannot get this TIFF image, which was exported from R to display in ArcMap 10.5. I would like to georeference the nodes, but the image only displays as all white.

Can someone help me with this?

I have tried removing the alpha layer in Gimp, exporting as black and white, resetting the gamma, changing the background display, using stretch instead of RGB, calculating statistics. The image renders a preview in ArcCatalog but not in ArcMap.


texts: Cell Statistics Error 000864: Overlay statistic: The input is not within the defined domain

texts:
I am trying to find the sum of a list of rasters using Cell Statistics in ArcPy. I am getting the following message: 'ERROR 000864: Overlay statistic: The input is not within the defined domain.'

I know this has something to do with not correctly coding the part where myrasters is input into the Cell Statistics function, but I can't figure out how to fix it. The function requires a list of rasters.

ESRI's sample code for this function is:

outCellStatistics = CellStatistics([inRaster01, inRaster02, inRaster03], "SUM", "NODATA")

Instead of writing out the whole list of rasters one by one for this function, I'm trying to get it to work with the myrasters list as an input.

Here's what I've got:

import arcpy
from arcpy import env
from arcpy.sa import *
import os

myrasters = arcpy.ListRasters('*.tif')
outraster = arcpy.gp.CellStatistics_sa([myrasters], "SUM", "DATA")
outraster.save('outworkspace\outname.tif')

UPDATE:
Switching to arcpy.gp.CellStatistics_sa(myrasters, "SUM", "DATA") as suggested by user2856 produces the same error. My list of rasters is structured as follows:

print(myrasters)
type(myrasters)
type(myrasters[0])

produces the output:

[u'raster1.tif', u'raster2.tif', u'raster3.tif', u'raster4.tif']
<type 'list'>
<type 'unicode'>


texts: How to anti-alias 3D projected contour lines, ArcScene is jaggy

texts:
I'd like to render some high quality contour lines in 3D. The image below was created in ArcScene. Note how ArcScene doesn't use any antialiasing, so the contour edges look jagged and nasty. Click on the image to see it at native resolution, it's only 736 pixels wide.

Does anyone know a solution to create high quality, antialiased GIS vectors in 3D. Can QGIS do it, or should I look into another solution like AutoCad or perhaps another 3D package? Ideally I'd rather stick to a GIS package than a non GIS 3D renderer as I don't have much familiarity with 3D file formats.

If anyone wants to access the 3D shapefile (750kb) it is here.


texts: Nearest neighbor analysis not showing results

texts:
I'm trying to do a nearest neighbour analysis on fish shoals distribution.
I go to  Vector > Analysis Tools > Nearest neighbour analysis. I choose the appropriate layer and when I hit ok the black and white wheel is on, and nothing happens; I think it IS doing the analysis because the blue bar is 100%; but just not showing anything.

Has anyone had this problem and know how to solve it?

I'm using: QGIS 1.8.0-Lisboa
GDAL 1.9.1
With Mac OS X 10.7.5


texts: Transform CRS of raster with PROJ4 string in pyqgis via QgsCoordinateTransform()

texts:
I try to reproject the crs of a raster to another common coordinate system, using this code:

crs = QgsCoordinateReferenceSystem()
crs_1 = crs.createFromProj4("+proj=tmerc +lat_0=0 +lon_0=38.03333333333 +k=1 +x_0=2400000 +y_0=-6511057.63 +ellps=krass +towgs84=23.57,-140.95,-79.8,0,0.35,0.79,-0.22 +units=m +no_defs")
crs_2 = crs.createFromProj4("+proj=utm +zone=37 +datum=WGS84 +units=m +no_defs")
QgsCoordinateTransform(crs_1,crs_2,QgsProject.instance())

But it gives me following error:

Traceback (most recent call last):
  File "[...]\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
TypeError: QgsCoordinateTransform(): arguments did not match any overloaded call:
  overload 1: too many arguments
  overload 2: argument 1 has unexpected type 'bool'
  overload 3: argument 1 has unexpected type 'bool'
  overload 4: argument 1 has unexpected type 'bool'
  overload 5: argument 1 has unexpected type 'bool'

What am I doing wrong? Using EPSG numbers works but is for me no option.

(Win 10, QGIS 3.1)


texts: gdal.Open error: 'NoneType' object has no attribute 'GetRasterBand'

texts:
I am running a python script which calibrates a hydrologic model, this script uses gdal from Osgeo. Following is the command I am giving python:

from osgeo import gdal
infile="G:\ARORA\vic_inputs\routing_files\new_basin\sutlej_beas_new_fraction.tif"
ds = gdal.Open(infile,gdal.GA_ReadOnly)
b1 = ds.GetRasterBand(band)

b1 gives the following error:

AttributeError                            Traceback (most recent call last)
<ipython-input-33-ba605a73e25c> in <module>()
----> 1 b1 = ds.GetRasterBand(band)
AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

The type(ds) returns NoneType.

What could be wrong here?


texts: Implementing OpenLayers3 basemap from tiles created using gdal2tiles

texts:
I recenty discovered the gdal2tiles.py tool which is working great with OpenLayers 2.x, but I wanted to use the new release of OL and I'm having some troubles. I have googled for a solution but i did not found a way to go.
Here is my code, inspired from the solutions offered here: https://groups.google.com/forum/#!searchin/ol3-dev/TMS/ol3-dev/JLn0lfHvI7w/q7C3O8_b0-UJ (still referring to the dev version of OL3, different it seems to the release).

var tileUrlFunction = function(tileCoord, pixelRatio, projection) {
  if (!tileCoord) { return ""; }

  var x = tileCoord[1];
  var y = tileCoord[2];
  var z = tileCoord[0];

  if (x < 0 || y < 0) { return "";}

  return 'http://localhost/map-tiles/'+ z.toString()+'/'+ x.toString() +'/'+y.toString() +'.png';
}

var layers = [
  new ol.layer.Tile({
    source: new ol.source.XYZ({
      tileUrlFunction: tileUrlFunction
    })
  })
];

var map = new ol.Map({
  target: 'map',
  layers: layers,
  renderer: 'canvas',
  view: new ol.View({
    center: [2900,-3900],
    zoom: 3
  })
});

I've tried several combinations of the x, y and z assignment without success. Tiles on the browser are completely wrong and or missing.


texts: Deleting Blank Spaces in Excel sheet before Importing to Table?

texts:
I have developed a model that imports an excel sheet and converts it into a table. It then joins values within this table to a particular feature class. I then use the calculate field tool to copy over the values in the table to the corresponding field in the Feature class. Only excel rows with values in each column of the excel table are copied over, if there is no value in the excel cell it does not copy over and leaves the original value that was in the feature class field (I'll attach image below). The issue is, many of the excel tables I've received are from multiple users and there are many instances of cells that are "blank" but there are spaces inserted in these areas in the sheet. The GIS see's these "blank" spaces with spaces as string values and will copy them over into the feature classes field it is joining. I need a front end process that will look at each excel spreadsheet first and delete out the spaces from cells that contain only spaces with no other actual value.

Is there any quick python script I can add into the model first to do this or to run after the excel is converted to a table?

This is a snippet of the model itself. The excel to table imports the particular sheet in the excel workbook and then joins it to the resulting feature class. The calc fields then copy over updated values from the spreadsheet to each of the feature classes corresponding field. At the end of the process, the imported fields from the excel sheet that were joined are then deleted. This is basically just a way to regularly update the feature classes as updates are made in the excel worksheets.

Here is the code I have so far for the calculator expression:

Expression:

testme(!ExcelTableFieldName!,!FeatureClassFieldName!)

Pre-logic script code:

def testme(x,y):
 if x is None:
  return y
 elif x == '':
  return y
 else:
  return x


texts:
I'm interested in wrapping long legend text (in ArcGIS 9.3.1) so that the text does not exceed the layout extent.  I would like to call the GetFormattedString Function within the UIButtonControl1_Click procedure.  Any pointers on how to get the legend text to wrap at 25 characters would be great.

Here is the code:

Public Sub UIButtonControl1_Click()

'Refresh legend

    Dim pMxDoc2 As IMxDocument
    Dim pPageLayout2 As IPageLayout
    Dim pGC2 As IGraphicsContainer
    Dim pElem2 As IElement
    Dim pMSF2 As IMapSurroundFrame
    Dim pMS2 As IMapSurround

    Set pMxDoc2 = ThisDocument
    Set pPageLayout2 = pMxDoc2.PageLayout
    Set pGC2 = pPageLayout2
    pGC2.Reset

    Set pElem2 = pGC2.Next
        Do Until pElem2 Is Nothing
        If TypeOf pElem2 Is IMapSurroundFrame Then
            Set pMSF2 = pElem2
            Set pMS2 = pMSF2.MapSurround

                Dim strFormattedString As String
                strFormattedString = GetFormattedString("TextContentblahblah", "Arial", 5, 12)

        pMS2.Refresh
        pMxDoc2.ActiveView.Refresh

        End If
        Set pElem2 = pGC2.Next

        Loop

End Sub

Private Function GetFormattedString(strTextContent As String, strFontName As _
String, sngFontSizeInPoints As Single, sngMaxWidthInPoints) As String

    Dim varWordArray As Variant
    Dim pMxApp As IMxApplication
    Dim pAppDisplay As IAppDisplay
    Dim pTransformation As ITransformation
    Dim pTextFont As IFontDisp
    Dim pTextSymbol As ITextSymbol
    Dim dblXSize As Double
    Dim dblYSize As Double
    Dim strGoodWidth As String
    Dim strFinalString As String
    Dim strTestString As String
    Dim i As Integer

    '------ Split the string into an array of words
    varWordArray = Split(strTextContent, " ", , vbTextCompare)

    Set pMxApp = Application 'QI
    '------ Get a pointer to the application display
    Set pAppDisplay = pMxApp.Display 'QI
    '------ Set up the Transformation
    Set pTransformation = pAppDisplay.DisplayTransformation
    '------ Set up the Font
    Set pTextFont = New StdFont
    pTextFont.Name = strFontName
    pTextFont.size = sngFontSizeInPoints
    '------ Set up the TextSymbol
    Set pTextSymbol = New TextSymbol
    pTextSymbol.font = pTextFont

    '------ Add each word into the test string and test for width
    pAppDisplay.StartDrawing pAppDisplay.hdc, 0
    For i = LBound(varWordArray) To UBound(varWordArray)
        If strTestString = "" Then
            strTestString = varWordArray(i)
        Else
            strTestString = strTestString & " " & varWordArray(i)
        End If
        '------ Get the TextSize
        pTextSymbol.GetTextSize pAppDisplay.hdc, pTransformation, strTestString, dblXSize, dblYSize
        '------ If the word added is < max width keep adding to the line, else make new  line by inserting vbcrlf
        If dblXSize < sngMaxWidthInPoints Then
            strGoodWidth = strTestString
        Else
            strTestString = varWordArray(i)
            strFinalString = strFinalString & vbCrLf & strGoodWidth
        End If

    Next i
    strFinalString = strFinalString & vbCrLf & strGoodWidth
    pAppDisplay.FinishDrawing
    'remove first vbcrlf and pass back
    GetFormattedString = Right$(strFinalString, Len(strFinalString) - 2)

End Function


texts:
I am dealing with long  and demanding analyses, which needs  often lots of data and time ( between few hours to couple of days) and it can happen a full stop/break on my workstations caused by diverse failures (data bugs, overflow, unexpected etc)

Is it possible to develop a  "Process Notification"  service, which can send an e-mail from ArcCatalog  to a chosen recipient about the actual status  , i.e.   Reporting  the status  in a brief form  : [1] Finished  [2] Break or Failure  [3] Estimated time left to be finished.

Our solution is not covering entirely the problem above .

For details:
[ http://geomatika.mapcat.net/COMMON/Projects/process_notifier_v0.41.zip ]


texts: Issues with getting the height from SRTM gdalwarp WGS84(EPSG:4326) to Mercator(EPSG:3857)

texts:
I am trying to use the SRTM .hgt file to obtain the height of a location. My issue is that I wish to reproject the SRTM file from WGS84 to Mercator but the height returned is not consistent after reprojection of the height file. From QGIS I selected a point where I believed that will cause such an issue.

The band value is different due to the fact that the coordinate reference system in QGIS can only read either WGS84 or Mercator at any point of time. This image is to ensure the  conversion of WGS84 coordinate to Mercator coordinate is correct when i am trying to retrieve the band value(height).

This image displays the issue I have been facing for quite sometime, as we can see that the coordinate have been converted to Mercator correctly, but the band value (height) returned is different from the WGS84, Mercator and Mercator with bilinear resampling.

What is the reason to cause the issue of different height value on different projection?

Command Used in Conversion:
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 N45E005.hgt MN45E005.tif
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 –r bilinear N45E005.hgt MBN45E005.tif

Downloadable Files to replicate test:
SRTM3 Files with created outputs


texts: Loading *.mrr files in QGIS or using GDAL

texts:
Does anyone know how to load the new proprietary MRR files using QGIS or GDAL?

I tried using Add Vector Layer but QGIS did not recognize this file.


texts: make Carto torque points change over time start from begin date and disappear at end date

texts:
First time use the torque function... I have variables: IDnumber,start_date, end_date, size, longitude, latitude.

Anyone know how to make carto torque points change depends on different dates?
I want to show the point on start dates, and make it disappear on its end date.

also I want to make each point size different depends on its size value.


texts: Viewshed analysis plugin on QGIS keeps resulting in errors

texts:
I've been trying to test the Viewshed analysis plugin on QGIS 2.10, but it keeps giving back error messages. I've tried to reduce the radius area and have also reprojected both my DEM and my point shapefile to UTM, in order to avoid a possible memory error mentioned in another question regarding this plugin (that the using of geographic coordinates might cause a memory error).

Still, I keep getting Python error messages when I try to run it. Does anyone know whats going on? Is this plugin "broken"?


texts: How to get a good resolution Image from google earth and use it in Arcmap?

texts: I need to make an orthophoto from google-earth. My orthophoto covers a large area and whenever I save the image in google-earth and put it in the Arcmap to rectify it , the image has very large resolution. I have found a software called "Google Satellite Maps Downloader" (which is not free). This software takes images for smaller areas from google-earth and stick them together, providing an ortho photo with better resolution. Unfortunately  my orthophoto is for a big area and the image which I get from the free version of this software  still has a bad quality. I was wondering if there is another way to get high quality satellite image from google-earth?


texts: Getting province and city name of given point

texts:
OpenStreetMap has many shapefiles for downloading, and with shp2pgsql we can import it to PostGIS

shp2pgsql -s SRID SHAPEFILE.shp SCHEMA.TABLE | psql -h HOST -d DATABASE -U USER

How can I get province and city name (or id) of a given point?
It is a little confusing for me! can someone explain more how can I do it? I just find some shape files here https://download.geofabrik.de/north-america/us.html and have no idea how to import just state/providence/city boundaries
Just tell me how to convert it to a minimal postgres table that contains type (state or city) , name of state/city and geometry of state/city
I know how can I use functions like ST_Contains to find city of given point


texts: Display Geojson data from PHP API (http://localhost/sample.php) on google map

texts:
I have PHP API (http://localhost/api.php) that gives data in GeoJSON format, how can I show that GeoJSON data on Google map ?

Dragging & Dropping GeoJSON data on map is not allowed.


texts:
We have one licence for Data-Interoperability server extension and one for Arcmap 10.6.

I want to share some of the tools that I created with some colleagues who do not have the extension. Basically, I want to run, share, and then allow them to connect to the server and run the tools with their own data.

I have a tool which has two parameters. Reader and writer. The reader is an SQL Server and the user must select the database which has to run through the whole flow. The output is a GDB with two polygon layers.

Problem 1:
From some reason the published tool doesn't have the list of SQL Databases.

Problem 2: The writer is static, even during publishing, and will not let me change it, no matter what.


texts: qgis atlas: represent only polygons that are within the main polygon

texts:
In the following image, the polygon without filling and the polygon with the red filling are two different layers. In my atlas I would like that only those red Polygons are represented that are inside the polygons with no filling.

The atlas is iterating over the polygons with no filling.


texts: ModelBuilder works in editor mode, not when opened?

texts:
I have made a pretty simple model. I can run it in the edit mode within ModelBuilder, and it works perfectly.

As you can see, I have set some parameters, so I do not need to run it in ModelBuilder. When I try to run the model then, it does not work.

The processing window says, that everything succeeded.

So what am I doing wrong with my model, since it can not work outside the edit mode?

I have read some others issues, and they should make a script. I want to avoid that.


texts:
We are hosting a cartodb installation ourselves which works pretty well but I can't seem to figure out how to get cartodb.js to call our installation instead of the one at cartodb.com.

Reading the source it's clear that cartodb.com is the default host so it should be possible to somehow specify another host.

The end goal is to build the maps programmatically using cartodb.js and the sql api (if that makes a difference).


texts:
I have a vector polygon layer which represents water bodies and I have a raster layer (250m pixel NDVI values). What I need is to mask all pixels in a raster layer that intersect my water mask and turn them into NoData.

Since I want to do it in ArcMap, there seems to be a bunch of tools like set null and con, but they work only with raster data. On the other hand, extract by mask does the opposite to what I need: it masks everything BUT the mask.

So my question is: is there a (best) way to apply a vector mask to a raster layer in ArcMap? And if no - how do I do it through rasterizing vector data, taking into account the resolution of the target layer?


texts:
I am working on a protected areas map. The label layer is from the mapping authority. The one I am using is point based.

Unfortunately some of the points carry the same name leading to duplicate names for mountain areas, glaciers and other. I am for several reasons bound to use this layer, so finding a different layer or making changes to the layer is not within the scope of this question.

Is there a way to remove duplicate labels when labels are within a certain distance of each other?

If so, can such a rule be enforced within one of the categories when rule-based labelling is used?


texts: vDatum returns error at Linux command line - fails to accept parameters with GUI

texts:
I am trying convert a GeoTiff using VDatum at the Linux command line.
My data reference the Washington State Plane South CRS with z values referencing EGM96.  I would like to convert these z values to NAVD88.
Here is my command.

/usr/lib/jvm/java-1.11.0-openjdk-amd64/bin/java -jar /opt/vdatum/vdatum.jar ihorz:NAD83_2011:spc:ft:4602 ivert:EGM1996:ft:height:GEOID12B ohorz:NAD83_2011:spc:ft:4602 overt:NAVD88:ft:height:GEOID12B -file:geotiff:geotiff,autosrcgeo:/home/gerry/geotiff/2019_1_308_405.tif;/home/gerry/geotiff/results/ region:3

I can confirm that I can run the vDatum GUI from the command line so  I know the paths to java and VDatum.jar are fine.  I know the paths to the geotiff are good.  This is the resulting error message:

NOAA's Vertical Datum Transformation - v4.1.2

System.file.separator= /
System.path.separator= :

cmdargs  :[Ljava.lang.String;@3930015a
12:37:13 : [vdatum_command]: java.lang.NumberFormatException
For input string: "/home/gerry/geotiff/2019_1_308_405.tif"

bash: /home/gerry/geotiff/results/: Is a directory

I do not understand why I would get a NumberFormatException on a file path string.
Thinking I would try to do the conversion manually I started the VDatum GUI.
Trouble with the GUI is that if I enter EGM96 for the input vertical datum, I get a pop-up claiming I need to reset the horizontal CRS.  If I reset the CRS to NAD83_2011 the vertical datum input reverts to NAVD88.  At the GUI vDatum will not accept the data vertical and horizontal values.
Thinking there might be a bug for the GUI Linux version of vDatum tried the input arguments at the online verson of vDatum. Online I can set the correct vertical and horizontal values for my input data but there is no geotiff conversion option online.  Now I am unsure if just the GUI version is malfunctioning, my command line is wrong, or if vDatum on Linux just does not work.
How do I use the Linux version of vDatum to convert Washington State Plane South, feet, EGM96 data to the NAVD88 vertical datum?


texts: Converting polygon to line feature in ArcGIS Desktop?

texts:
I received a data file containing city sidewalks as a polygon feature. I am running ArcGIS 10.3.1, with a Spatial Analysis extension. For the purpose of my data analysis, I need these polygon features to be single line features.

I used the "Polygon to Line" tool, which produced a line feature of the perimeter of each polygon (not what I'm looking for). From this output, I used the "Collapse Dual Lines to Centerline" tool, with a maximum width of 1 foot. This basically produced the same output as the input, with the exception of two ouliers (not centerlines).

Can anyone provide a solution or suggest a workflow?


texts: Coordinate system of a MODIS file to be introduced in GDAL for transformation

texts:
I am working on MODIS data using GDAL. When I extract coordinate system information of the files I got this:
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Unknown datum based upon the custom spheroid",
        DATUM["Not specified (based on custom spheroid)",
            SPHEROID["Custom spheroid",6371007.181,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (2223901.039533000000000,0.000000000000000)
Pixel Size = (926.625433139166600,-926.625433139166720)

I have to transform the coordinate system of my other files into this coordinate system. But the problem is that I don't know how to introduce this coordinate system to my program.
I can't find EPSG or Proj4 for it. Any solution or alternative?


texts: Sentinel-1 GRD and Sentinel-2 L1C on AWS bucket policy

texts:
Recently I started to work with Sentinel-1 data hosted on AWS S3 bucket 'sentinel-s1-l1c'. The bucket is requester pays and requires AWS signature version 4. I need to access the files via HTTP GET requests, so I generate authorization headers myself as described in AWS docs.

When the signature expires there's an XML exception, saying that MaxAllowedSkewMilliseconds equals 900 seconds. I estimate the expiration time as date in X-Amz-Date header (time of the request) + 900 seconds. But this method seems incorrect, since I have 403 Forbidden earlier than expected.

What is the correct method to calculate signature expiration time?
Should I consider SignatureAge in bucket policy? If so, then how can I get the policy? I tried to get the policy using boto3 and awscli but no success...

Since Sentinel-2 also goes requester pays, I wonder what the policy will be and how to determine the signature age when migrating code to consider Sentinel-2 bucket as requester pays?


texts: Dictionary (Error) Error in map(ID=4): Date: Parameter 'value' is required. Error in Google Earth Engine

texts:
I am using the savizky Golay method to smooth NDVI data. My data is Landsat NDVI generated by the Maximum Value Composite method. My code reference from https://code.earthengine.google.com/e9f4f2bb84f0bdd26ecfbea47a71885f, and my code is at https://code.earthengine.google.com/9c06176826fca41c06eba08a64ae9bf9.  My code runs with an error:

Dictionary (error) error in map (id = 4): Date: parameter 'value' is
required.

I can't find a way to solve it. I hope someone can help me.


texts: st_transform producing spurious ribbons across region

texts:
I am trying to use st_transform on PERZONES shapefile available here . In its native projection, it looks like this:

However, when I run the command:

tansformedPF<-st_transform(PF,"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
##This should be the MODIS sinusoidal projection

I get:

I thought maybe st_transform_proj would fix the problem but it did not. I also tried to plot the transformed file in Earth Engine just in case it was a problem with R's plotting and not the raster itself, but EE also shows the bands. I've used st_transform with the same destination crs with other shapefiles with no problem.

I also tried st_wrap_dateline(PF, options = c("WRAPDATELINE=YES") and st_wrap_dateline(PF, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")  before passing to st_transform, but no luck.   The native proj4string is "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs".

Any workarounds to reproject this file without the bands?  I am open to another reproject function, but st_transform seems to be the one to use.


texts:
I am using ArcMap 10.6.1.
I have a GeoJSON file and I wanted to load into my map. I click on the add data option and go to the location where geojson file is stored. However, I could not see the geojson file located on the folder.
Am I doing this in the right way?


texts:
I'm trying to create a template map following this tutorial from the cartoDB's maps API.

My template.json file

{
  version: "0.0.1",
  name: "template_1",
  // embedded authorization certificate
  auth: {
   method: "open"
  },
  placeholders: {
      color: {
        type:"css_color",
        default:"blue"
      },
      del_s: {
        type:"sql_literal",
        default: "19420210"
      },
      del_n: {
        type:"number",
        default: 19420210
      }
  },
  layergroup: {
    "version": "1.0.1",
    "layers": [{
      "type": "cartodb",
      "options": {
        "cartocss_version": "2.1.1",
        "cartocss": "#layer { polygon-fill: <%= color %>; }",
        "sql": "select * from table1 WHERE del = <%= del_s %>",
        "interactivity": ["cartodb_id","del","gk"]
      }
    },
    {
      "type": "cartodb",
      "options": {
        "cartocss_version": "2.1.1",
        "cartocss": "#layer { polygon-fill: red; }",
        "sql": "SELECT * from table2 WHERE del = <%= del_n %>",
        "interactivity": [ "cartodb_id","del"]
      }
    }]
  }
}

When I try to make the request to create the template..

curl -X POST \
   -H 'Content-Type: application/json' \
   -d @template.json \
   'https://docs.cartodb.com/api/v1/map/named?api_key=568kaka8bc89'

I can not connect, I get an error:

curl: (7) Failed connect to docs.cartodb.com:443; Connection timed out

anyone knows what I'm doing wrong?

UPDATE

Try with my username, and get an "internal server error"

curl -X POST \
   -H 'Content-Type: application/json' \
   -d @template.json \
   'https://myusername.cartodb.com/api/v1/map/named?api_key=myapikey'


texts: Creating LineString and length of LineString from multiple latlon points for each row in a Pandas Dataframe

texts:
So I want to create a LineString and get the length of the LineString from two sets of LatLon coordinates (points)  in each row in a Pandas Dataframe, and store in a new column

How can I do this?


texts: Dissolve polygons based on number of polygons

texts:
How can I dissolve group of polygons into many polygons based on specific number?

For example I have 100 points which are beside each other I need to dissolve each nearest 5 polygons into 1.


texts: Annotation group for data frame grid in ArcGIS

texts:
I am working in a project which I should label two types of grids (one is projected and the other is geographic) in the same data frame (layout view).

Often there is an overlap of these coordinate labels. I usually have to convert one of them to graphic (Data Frame Properties > Grids > Select Graticule grid > Convert to Graphics) to manually adjust them.

When I am working in the next geological sheet (I am using the same layout template for the different sheets) I need to manually delete the grids features converted to graphics to start again.

I would like to know if this stored graphical grid is located somewhere. For example, when I convert some label to annotation from a shapefile I should find it in the Annotation Groups (Data Frame Properties) and erase it. Is there a similiar option when you convert a grid to graphics? There is only the Default Group Name in the Annotation Groups.

PS: My ArcMap version is 10.3.


texts: I'm working on a project of water supply which includes a layer representing the pipes. the operation that I want to do is generating a point layer where a pipe touches an other and then split it at this location. I found out that st_touches in PostGIS does it but I'm wondering if I can do it on the same layer because in the doc they were talking about two different layers or geomtries. How can I achieve this?


texts: Why wont QGIS let me export selected features only?

texts: I am trying to export a specific set of soils from a state wide soil survey. However, when I try use the export selected features it wont let me use the feature?


texts: Most Efficient method of Storing Point Data in PostGres/PostGIS to later re-create a track

texts:
I am trying to determine the most efficient way to store point data in PostGIS to later re-create tracks

  Criteria:

  1 second updates

  100- 200 tracks active at any one time

  1 hour max track time so approx 5400 points per track

  Track data extracted as a report, not real time display, but the reports would be run frequently, many times daily

I can see two possibilities at this stage:

    a) Store the point data as a table with a field that will group them to a particular track. So three fields used, gid, trackid, lat/lon, resulting in a lot of rows

    or:

    b) Store each track as one row in a table with a field that has the point data appended to it at each update

My Questions are:

  Which is a more efficient way to store the point data, Numeric Fields (two columns), char fields, or GIS POINT/LINESTRING fields
  Which would be the most efficient way of the two options above, given server/database load to store and recreate the track data. Option a) would require queries to build each reported track whereas option b) would only need to look at one row per track.
  Is there a better way than either a) or b) to approach this requirement.

Thanks
Mark


texts: Rotating markers for points in Tilemill along direction of a way

texts:
I have designed a simple marker to represent pedestrian crossings, which in OpenStreetMap data are represented as nodes (points) along the way (line). I want to rotate the marker to follow the line it sits on. But I can't see a way to do that.

Is there a way?


texts: OpenLayers 2.13.1 - edit and export GeoJSON property values

texts:
I have created an application with OpenLayers 4.1.0 that imports GeoJSON vector data, allows selecting of polygons, changes "processed" property value of selected polygon from 0 to 1 and prepares edited data for export.

Code is here:
https://github.com/mleontenko/openlayers-edit-properties/blob/master/map.html

Now, I'm trying to rewrite that application in OL 2.13.1:
https://github.com/mleontenko/openlayers2-edit-properties/blob/master/map.html

This is the part i can't figure out how to rewrite in OL2:

//change "processed" property value from 0 to 1
select.getFeatures().item(0).set('processed', 1);

//store geojson to variable for export
var writer = new ol.format.GeoJSON();
var geojsonStr = writer.writeFeatures(vectorSource.getFeatures());enter code here

What it needs to do is:

-Change selected polygon's "processed" property from 0 to 1

-store GeoJson to variable as string so i can export it later


texts: ST_Area(geometry) with ogrinfo estimates a stable area per grid cell

texts:
I use ST_area(geometry) in order to calculate the area of each grid cell of a vector grid with LCC projection and 6x6km2 resolution. The CRS of the vector is:
+proj=lcc +lat_1=45 +lat_2=22 +lat_0=38.7 +lon_0=20 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs

I do the following:

ogrinfo -dialect SQLite -sql "SELECT id,ST_Area(geometry) as AREA FROM grid" grid.shp
ogrinfo grid.shp -sql "ALTER TABLE grid ADD COLUMN AREA_grid float"
ogrinfo grid.shp -dialect SQLite -sql "UPDATE grid SET AREA_grid =ST_Area(geometry)"

The result is that the area_grid takes the value of 36000000m2 for all grid cells. This is not correct since the real area of grid cell slightly varies in different coordinates.

Do I have to define in a different way the area of the grid cell?


texts: Reading geometry that starts with "0xB508000..." and may be WKB

texts:
Someone has given me an Excel file that seems to have been exported from ArcGIS. All columns are fine, except for the geom field. Every row in that field starts with 0xB508000... and some number of zeroes.
ST_GeomFromWKB and similar PostGIS functions all yield errors.
Here's a sample. It should yield a parcel somewhere on Catalina Island.
0xB5080000010405000000C01E8567E7A55841C0CCCC845D26384150B81EF5DDA55841007B14D67A263841D0A37094D8A55841C047E1065F263841608FC20BE2A55841000000D441263841C01E8567E7A55841C0CCCC845D26384101000000020000000001000000FFFFFFFF0000000003


texts: "Parse error at line 1 column 1" when using QConsolidate3 to share QGIS projects

texts:
I am trying to save a project and layers as a zip file to share with colleagues but I am receiving an error message when I do so:

  Parse error at line 1 column 1

I am using QGIS desktop 3.2.3 and QConsolidate3.


texts:
I tried Google, but the questions isn't that easy to put in Google, so I haven't found anything about this.
I am making a database with different locations (latitude and longitude).
For this I have two GPS systems, a differential and a non-differential. My question is, can the positions measured by the different systems be compared?
As I understand, the accuracy is only a statistical value.  So if the accuracy tells you 5 cm, you know that the systems are measuring in this radius. But how do you know that the measurements are correct to begin with?

When I complete the database, I can compare further measurements made with it to check the performance of the GPS.  Is it then better to make the database with the differential, or make two different ones because I can't compare the two systems with each other?


texts:
I am new to postgis and I need to find intersection of three polygons.
In my table named "mytable" contains large no of polygon coordinates and grouped into three categories by adding a field called "filename".I labelled the category names as part1, part2, part3 in filename . Now I need to find all polygon which spacially intersect in three categories.
I write my code as:

SELECT ST_Intersection(a.geom, b.geom) as intersect_ab
FROM mytable a
INNER JOIN mytable b
ON ST_Intersects(a.geom,b.geom)
INNER JOIN mytable c
ON ST_Intersects(a.geom,c.geom)
WHERE ST_Overlaps(a.geom, b.geom)
AND ST_Overlaps(b.geom, c.geom)
AND ST_Overlaps(a.geom, c.geom)
AND ST_isvalid(a.geom)='t'
AND ST_isvalid(b.geom)='t'
AND ST_isvalid(c.geom)='t'
AND a.filename = 'part1'
AND b.filename = 'part2'
AND c.filename = 'part3');

But it returns wrong answer


texts: Defining a coordinate system in WKT or PROJ format that has an Affine transformaiton and bounds

texts:
I need to define a local coordinate system via OGC WKT or PROJ form.

I know they have implemented Affine transformations in Proj 5 but I am unsure of how to implement it. The base projection is from EPSG 28354 but we have file in the local coordinate systems with the below Affine transformation and Bounds

Name    SMG
Projection  Transverse Mercator, (also known as Gauss-Kruger)
Datum   GDA 94 Australia GRS 80
units   meter
Origin Longitude    Central Meridian
Origin Latitude used    central Parallel
Scale Factor    0.9996
False Easting   500000
False Northing  10000000
Affine Unit meter
Affine A    0.88866511
Affine B    -0.459520111
Affine C    3099033.308
Affine D    0.459520299
Affine E    0.888664756
Affine F    -6949726.605
Bounds  MinX    -10000
Bounds  MinY    -10000
Bounds  MaxX    500000
Bounds  MaxY    500000

I tried defining it with FITTED_CS but QGIS (our GIS software does not recognise that). Here is what I have tried:

FITTED_CS["test",
 PARAM_MT["Affine",
  PARAMETER["num_row",3],
  PARAMETER["num_col",3],
  PARAMETER["elt_0_0", 0.88866511],
  PARAMETER["elt_0_1", -0.459520111],
  PARAMETER["elt_0_2", 3099033.308],
  PARAMETER["elt_1_0", 0.459520299],
  PARAMETER["elt_1_1", 0.888664756],
  PARAMETER["elt_1_2", -6949726.605]],
 PROJCRS["MGA",
  BASEGEODCRS["GDA94",
    DATUM["Geocentric Datum of Australia 1994",
      ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1.0]]]],
  CONVERSION["Map Grid of Australia zone 54",
    METHOD["Transverse Mercator",ID["EPSG",9807]],
    PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.01745329252]],
    PARAMETER["Longitude of natural origin",141,ANGLEUNIT["degree",0.01745329252]],
    PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1.0]],
    PARAMETER["False easting",500000,LENGTHUNIT["metre",1.0]],
    PARAMETER["False northing",10000000,LENGTHUNIT["metre",1.0]]],
  CS[cartesian,2],
    AXIS["easting (E)",east,ORDER[1]],
    AXIS["northing (N)",north,ORDER[2]],
    LENGTHUNIT["metre",1.0]]]


texts: Why are some QGIS tools only available via the Processing Toolbox?

texts:
This is a followup to Making QGIS toolbox tool available directly from toolbar? which I posted several days ago.

That question was more of a "how".  This one is more of a "why".

Why are some QGIS tools only accessible via the Processing Toolbox, and not available directly via the normal menu?

For example, at 2.18.16 I constantly use the Frequency Analysis tool.  But every time I use that tool I am required to click Processing > Toolbox > QGIS geoalgorithms > Vector table tools > Frequency Analysis.

Whew!  It would be a lot easier if that tool was directly available on the main menu rather than hidden deep within the Processing Toolbox menu structure.

I'm not talking about plugins, but the tools that come with a default QGIS install.


texts:
I am trying to put the Solar Radiation data into QGIS with my current projection of EPSG:27700 but even when I change the projection, the data is well off the mark. The data comes from here: http://re.jrc.ec.europa.eu/pvgis/download/solar_radiation_classic_laea_download.html

The website states the data is in the map projection: Lambert Azimuthal Equal Area, ellipsoid WGS84


texts: Changing properties of QWidget using python init function with QGIS?

texts:
I have a form created with the Drag and Drop Designer in QGIS, and it contains a Relation Reference widget.

Is there any way to change the available XML configuration from within a python init function for the QGIS form?

<edittype widgetv2type="RelationReference" name="tree_uid">
      <widgetv2config OrderByValue="0" fieldEditable="1" AllowAddFeatures="1" ShowForm="1" Relation="test20160803125723549_tree_uid_taxa_tree_detail20160812111449449_UKEY" ReadOnly="1" constraint="" MapIdentification="0" labelOnTop="1" constraintDescription="" notNull="0" AllowNULL="0"/>


texts: Why is ArcGIS crashing when I try to load .adf raster data?

texts: I am pretty new to GIS, so forgive me if this is a stupid question.  I downloaded a bunch of land cover files from USGS seamless warehouse of Florida.  Every time I try to 1) view the .adf file (which I am assuming is the raster data image) in arcCatalog or 2) add the .adf file in arcMap the program crashes.  Any reason for this?  I am running the trial version of ArcMap 10 on my home PC in order to avoid using the slow ones at my universities lab, so I figured maybe I am missing some kind of extension or something.  I am not really sure, I do not have much experience working with raster data so far.


texts:
I have to mosaic several pairs of images (A and B).
"A", "B" images has the same format (tif), projection, cell size, extent.
Cutline in shape format is available if its needed.
Nodata values are not present.

Im looking for a scriptable (preferable gdal) solution to do this, but nothing useful found yet.


texts: Spatial index info does not appear in ALL_SDO_INDEX_INFO

texts:
I have a table with this definition:

CREATE INDEX "MYSCHEMA"."MYTABLE_SX" ON "MYSCHEMA"."MYTABLE" ("GEOLOC")
INDEXTYPE IS "MDSYS"."SPATIAL_INDEX"
PARALLEL 10 ;

Geoloc It is a spatial column, but I don't see it in the ALL_SDO_INDEX_INFO table

The

PARALLEL 10 ;

Has something to do? or maybe my column is missing something?


texts: QGIS API 3.2 and PyQT5 in Ubuntu Fails: Native Qt signal is not callable

texts:
I'm trying to migrate from QGIS API 2.18 to 3.2 in Ubuntu18. Several changes are made, including the update from PyQT4 to PyQT5. This is the very beginning of my code:

from qgis.core import QgsApplication, QgsProject
from qgis.gui import QgsMapCanvas,  QgsLayerTreeMapCanvasBridge
from PyQt5.QtCore import QFileInfo
from PyQt5.QtXml import QDomDocument

reportDir = ''
project_path = reportDir+'/miQgsVegProject2.qgs'
template_path = reportDir+'/miQptVegProject2.qpt'

QgsApplication.setPrefixPath("/usr/share/qgis", True)
app = QgsApplication([], False)
app.initQgis()

canvas = QgsMapCanvas()
project = QgsProject.instance()

project.readProject(QFileInfo(project_path))
layers=QgsProject.instance().mapLayers()
bridge = QgsLayerTreeMapCanvasBridge(project.layerTreeRoot(), canvas)
bridge.setCanvasLayers()
template_file = file(template_path)
template_content = template_file.read()
template_file.close()
document = QDomDocument()
document.setContent(template_content)
composition = QgsComposition(canvas.mapSettings())
composition = QgsComposition(canvas.mapRenderer())
composition.setPlotStyle(QgsComposition.Print)

composition.loadFromTemplate(document, {})
template_file = file(template_path)
template_content = template_file.read()
template_file.close()
document = QDomDocument()
document.setContent(template_content)
composition = QgsComposition(canvas.mapSettings())
composition = QgsComposition(canvas.mapRenderer())
composition.setPlotStyle(QgsComposition.Print)

composition.loadFromTemplate(document, {})

map_item = composition.getComposerItemById('mapa')
map_item.setMapCanvas(canvas)
map_item.zoomToExtent(canvas.extent())
#map_item.zoomToExtent(miExtent)

composition.refreshItems()
composition.setPlotStyle(QgsComposition.Preview)
composition.setPrintAsRaster(True)
composition.exportAsPDF(pdfOutput)
QgsProject.instance().clear()

app.exitQgis()

This tiny code results in:

  project.readProject(QFileInfo(project_path)) TypeError: native Qt
  signal is not callable

At PyQT4 this works but I can't find a clue to migrate correctly with PyQT5.


texts: How to stop transfer of large raster into arcsde geodatabase

texts:
I have just realised that I will not have enough space on my harddrive to accomodate a raster which is being copied into an arcsde geodatadatabase. I am afraid, that if I stop the process in ArcCatalog, I will corrupt the raster in the geodatabase and not be able to delete it...or even worse corrupt the whole database.

Can anyone advise me on how to stop the copying process and ensure that the sde Geodatabase will not hold any remnants.

Arcsde9.3 / ArcCatalog 10.0 / MSSQL


texts:
I have a shapefile with the attribute "yield". I need to know the math difference between yield (id1) and yield (id2)... then yield (id2) and yield (id3) and so on... I hope you understand what I mean.

In Excel it is easy to handle. But I do not want to export and import the files.


texts: Visualizing speed of movement on static (= printable) map?

texts:
Imagine you're tracking, for example, vehicle or animal movement and receive regular GPS position updates. How would you visualize speed of movement on a static (= printable) map? I've been coloring the points in red-yellow-green where red would be "slow", but I think there has to be a better/more intuitive way to visualize speed measurements.

Another requirement would be that it should be possible to track multiple moving objects at once. Their paths can cross or run alongside each other (on the same roads).


texts: How to find the frequency of pixel categories inside a geometry (Google Earth Engine)

texts:
The pixels of a certain Image are classified into 4 categories. I want to count how many pixels within a certain geometry fall in each category.

I thought my code below was doing pretty much that. However, there are two things about the results that I find very odd.

Thing #1: the counts have decimal places. They should be integers, no? They are, after all, counts.
Thing #2: I have two geometries with the same area. The number of pixels inside them should be, at least to an approximation, the same. That, however, is not the case. The total number of "pixels" in one geometry is more than twice the number of pixels in the other geometry.

These two facts lead me to believe that I am not achieving what I am trying to achieve. I'm stuck.

Am I achieving what I am trying to achieve?

Here's my code. You can run it in Google Earth Engine following this link. At the end of the code, I print the frequencies of each category in each geometry. In the first, the total frequency is 24 + 41616.96 +3019.09 + 26091.79 = 70751.85. In the second, we get 25.20 + 1706.94 + 1715.46 + 28427.35 = 31874.97.

// imports
var jrc2003 = ee.Image("JRC/GSW1_0/YearlyHistory/19"),
    pt1 = /* color: #d63000 */ee.Geometry.Point([-64.64, -9.27]),
    pt2 = /* color: #98ff00 */ee.Geometry.Point([-51.77, -3.13]);

// select the image band
var water = jrc2003.select(['waterClass'])

// create buffers
var bf1 = pt1.buffer(10 * 1000);
var bf2 = pt2.buffer(10 * 1000);

// areas of the buffers I just created -- they are the same
print('Area of bf1 in km2:', bf1.area().divide(1000 * 1000));
print('Area of bf2 in km2:', bf2.area().divide(1000 * 1000));

// this is what I think I should do to get the count of pixels,
// within each geometry, in each category
var dic1 = water.reduceRegion({
  reducer  : ee.Reducer.frequencyHistogram(),
  geometry : bf1,
  scale:30,
  maxPixels: 1e9});

var dic2 = water.reduceRegion({
  reducer  : ee.Reducer.frequencyHistogram(),
  geometry : bf2,
  scale:30,
  maxPixels: 1e9});

// print histograms -- results look odd to me
print('Hist 1', dic1);
{
  "waterClass": {
    "0": 24,
    "1": 41616.96470588233,
    "2": 3019.090196078431,
    "3": 26091.796078431376
  }
}

print('Hist 2', dic2);
{
  "waterClass": {
    "0": 25.20392156862745,
    "1": 1706.9450980392155,
    "2": 1715.4666666666667,
    "3": 28427.356862745102
  }
}


texts:
I have many SHP tables organised by year, for example 2015_xx_data.shp that I have loaded into a MS SQL Server database using OGR2OGR. I have also loaded them into MySQL for a comparison. There are no problems with the MySQL database but I am unable to get any features from the MSSQL database using OGRINFO.

I am using SQL Server 2014 Management Studio.

I understand that the table names really should not start with a numeric and would need to be qualified with quotes or [] in this case.
But the GDAL MSSQL driver does not seem to work quite right here.

QRY 1: If I use this syntax I get the Feature Count 0 and the correct SRS being reported, the fields being reported but no data is output:

ogrinfo -al "MSSQL:server=xxx\xxx;database=spatialtest;trusted_connection=yes;tables=dbo.2015_xx_data"

INFO: Open of `MSSQL:server=xxx\xxx;database=spatialtest;trusted_connection=yes;tables=dbo.2015_xx_data'
      using driver `MSSQLSpatial' successful.

Layer name: 2015_xx_data
Geometry: Line String
Feature Count: 0  <<<<<< !!!!
Layer SRS WKT:
PROJCS[...etc...]
FID Column = ogr_fid
isovalue: Real (6.2)

QRY 2: If I use this syntax I get Feature Count 3095 (correct) but the the SRS is (unknown) and no features are dumped because dbo.geometry_columns has "2015_xx_data" as the entry, not "[2015_xx_data]".

ogrinfo -al "MSSQL:server=xxx\xxx;database=spatialtest;trusted_connection=yes;tables=dbo.[2015_xx_data]"

ERROR 1: No column definitions found for table '[2015_xx_data]', layer not usable. <<<<<< !!!!
INFO: Open of `MSSQL:server=xxx\xxx;database=spatialtest;trusted_connection=yes;tables=dbo.[2015_xx_data]'
      using driver `MSSQLSpatial' successful.

Layer name: [2015_xx_data]
Geometry: Unknown (any) <<<<<< !!!! should be Line String
Feature Count: 3095
Layer SRS WKT:
(unknown)

So I manually added a new row into geometry_columns with the f_table_name "[2015_xx_data]". QRY 2 now returns Feature Count 3095, the correct SRS but no data is output. ERROR 1 shows the problem.

ogrinfo -al "MSSQL:server=xxx\xxx;database=spatialtest;trusted_connection=yes;tables=dbo.[2015_xx_data]"

ERROR 1: No column definitions found for table '[2015_xx_data]', layer not usable. <<<<<< !!!!
INFO: Open of `MSSQL:server=xxx\xxxx;database=spatialtest;trusted_connection=yes;tables=dbo.[2015_xx_data]'
      using driver `MSSQLSpatial' successful.
Layer name: [2015_xx_data]
Geometry: Unknown (any) <<<<<< !!!! should be Line String
Feature Count: 3095
Layer SRS WKT:
PROJCS[...etc...]
No data is output <<<<<< !!!!*

I have tried many permutations on this, no luck. But MySQL works flawlessly.
I think the MSSQL driver isn't behaving as it should here, possibly a bug?
Any advice on how to overcome this?

EXTRA INFORMATION

I have tried this with GDAL (ogrinfo.exe) 1.8, 1.11 and 2.0 - all same behaviour.

Quantum GIS Pisa 2.10.1 works fine and allows me to add the layers and they draw as expected. This means that QGIS or the GDAL library shipped with it is doing something different - not sure what.


texts: Contour extraction from raster does not finish in Qgis. How to find out if a raster is damaged

texts:
Execution of "extraction" - "contours" in Qgis is impossible.
The process never finishes, I stopped it after 35 min.

The contour extraction works for all other raster files.
I have tried with copies and converted the format but nothing works.

I can find the contour folder and the shapefile of the created contour but when imported to Qgis they don't show anything.

Is there a way to find out if the raster is damaged, or what cause there might be to this problem?


texts:
I am classifying Landsat 8 and Sentinel 1 in Google Earth Engine but the console tells me I do not have the required bands for classification even after assigning all the bands.
Attached a link to my workflow (JavaScript for GEE):
var season1 = Landsat.filterDate('2018-10-01', '2018-12-31')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season1:', season1);

var season2 = Landsat.filterDate('2019-01-01', '2019-03-31')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season2:', season2);
var season3 = Landsat.filterDate('2019-04-01', '2019-06-30')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season3:', season3);
var season4 = Landsat.filterDate('2019-07-01', '2019-09-30')
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 25))
                  .filterBounds(mis)
                  .map(maskL8sr)
                  .median().clip(mis);
print('Landsat_season4:', season4);

// Define False color visualization
var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B5', 'B4', 'B3'],
  gamma: 1.7
};

///Visualize the FCC of each season Landsat images
Map.addLayer(season1, rgbVis, 'Landsat MSI - season1',false);
Map.addLayer(season2, rgbVis, 'Landsat MSI - season2',false);
Map.addLayer(season3, rgbVis, 'Landsat MSI - season3',false);
Map.addLayer(season4, rgbVis, 'Landsat MSI - season4',true);

// Below four lines select season wise Sentinel 1 RADAR data again median and clipped to MIS
var vh_s1 =  S1.filterBounds(mis).filterDate('2018-10-01', '2018-12-31').median().clip(mis).select('VH').rename(['vh_s1']);
var vh_s2 =  S1.filterBounds(mis).filterDate('2019-01-01', '2019-03-31').median().clip(mis).select('VH').rename(['vh_s2']);
var vh_s3 =  S1.filterBounds(mis).filterDate('2019-04-01', '2019-06-30').median().clip(mis).select('VH').rename(['vh_s3']);
var vh_s4 =  S1.filterBounds(mis).filterDate('2019-07-01', '2019-09-30').median().clip(mis).select('VH').rename(['vh_s4']);

// Here we combine all the seasonal Landsat and Sentinel 1 to a single image collection called
// 'dataset'
var dataset = season1.addBands(season2).addBands(season3).addBands(season4).addBands(vh_s1).addBands(vh_s2).addBands(vh_s3).addBands(vh_s4);
print('dataset:', dataset);

//Below codes split the training data into Training (70%) and Validation (30%) - Random
var trdata_table = field.randomColumn({seed: 1});
var training = trdata_table.filter(ee.Filter.lt('random', 0.70));
var validation = trdata_table.filter(ee.Filter.gte('random', 0.70));

// Create a new feature collection (vector) with all the bands in 'dataset' sampled to training points
var training = dataset.sampleRegions({
  collection: training,
  properties: ['class_id'],
  scale: 20,
  tileScale: 8
});
print('training:', training);
// Create a new feature collection (vector) with all the bands in 'dataset' sampled to validation points
var validation = dataset.sampleRegions({
  collection: validation,
  properties: ['class_id'],
  scale: 20,
  tileScale: 8
});

//Define bands require for training and classification. So Bands 1,2,3,4,5,6,7,9,10,11 for 4 seasonal
// images and Sentinel 1 (1 band) representing 4 seasons
var bands = ['B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B9', 'B10', 'B11', 'B1_1', 'B2_1', 'B3_1', 'B4_1', 'B5_1', 'B6_1', 'B7_1', 'B8_1', 'B9_1', 'B10_1', 'B11_1', 'B1_2', 'B2_2', 'B3_2', 'B4_2', 'B5_2', 'B6_2', 'B7_2', 'B8_2', 'B9_2', 'B10_2', 'B11_2', 'B1_3', 'B2_3', 'B3_3', 'B4_3', 'B5_3', 'B6_3', 'B7_3', 'B8_3', 'B9_3', 'B10_3', 'B11_3', 'vh_s1', 'vh_s2', 'vh_s3', 'vh_s4'];
print('bands:', bands);


texts: How do I get ArcGIS tiles to display a single layer when multiple layers are available?

texts:
I am using Bing Maps v7 to display layers from an ArcGIS server. The server displays multiple layers as one, but I would like to only display a single layer.

Suppose a map service has the layers: Topographic Info (0), Places Info (5), and Citations(10)

I write code to get a tile like so:

var source = new Microsoft.Maps.TileSource({
        uriConstructor: function( tile ) {
            return baseUrl + "/tile/" + tile.levelOfDetail
                        + "/" + tile.y + "/" + tile.x;
        }
    });

Which would procdure a URL like this:
http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/12/1528/1218

Is it possible to have the generated tile only display information regarding a single layer, such as Topographic Info (0)?

I tried variations of things like ?layers=0 in the query string, but nothing seems to work.


texts: Maintaining Order and Repetition of Cell Names using PostGIS?

texts:
I am trying to get cells names (in a sequence) by way trajectory is passing through them. For example in the photo shown below:

Blue markers represent a grid and red line represent a trajectory passing through them. (Its a 5x5 grid).

The result should be{C22, C21, C31, C30, C40, C41} where C21 etc are the cell names by row & column number. All it means is trajectory is starting from cell_names C22 to C21 to C31 and so on.
I am using following query to retrieve this result:

select ce.cell_names
from cells ce, traj tr
where st_intersects(tr.traj_path, ce.coordinates) and tr.traj_id = 207 and ce.grid_id = 7776

Where traj_path is LINESTRING and coordinates is POLYGON.

The problem of keeping cell order was solved by the answer of one previous contributor and it also works through following query:

SELECT tr.traj_id, ce.cell_id
from cells_new_york AS ce, real_traj AS tr,
ST_LineLocatePoint(tr.traj_path, ST_Centroid(ce.coordinates)) as saim
where st_intersects(tr.traj_path, ce.coordinates) and ce.grid_id = 9 and tr.traj_id = 572
order by saim

But this misses the revisited polygons. As in the figure below C22 and C21 are visited twice. But they are not captured by the query.

Can you guide me towards a relevant PostGIS function?


texts:
I access this WMS (link) which is returned in only .png/.jpg/.gif in a pre-defined colour-scale. I want to change those colours and am asking myself if I can set up a custom style for the rasters. I'm now assuming the answer would be a bald 'no, it's a raster'. But I'm looking for maybe a work-around? I can also access the point-features which it is based on (link). Can I use some sort of styling or filters to colour the cells according their elevation ('on-the-fly')?

For more clarity I have the WMS running on here. The web page is in Dutch and the layer you should look for is 'Lodingen' (enabled by default).


texts: How to programmatically set Selection-> Selection Options [Scale selection symbols when a reference scale is set] checkbox checked or unchecked

texts:
I have seen the unanswered  Scale Selection Symbols with Reference Scale using ArcGIS Engine? at the time of posting this question, and my question is different as i need to set this "Scale selection symbols when a reference scale is set" setting in ArcMap version 10.2.1

I have seen couple of Interfaces related to Selection but none of those have any option to set this setting specifically.
below are the few links i have seen:

IMxApplication.SelectionEnvironment
ISelectionEnvironment
ISelectionEnvironmentStorage
ISelectionEnvironmentThreshold


texts: How can I show pictures from map tips in qgis2web exported maps

texts:
I have a point layer which pictures. The pictures are stored in my hard drive. In QGIS the map tips and feature attribute form show picture of point perfectly as shown in these 2 prin shots.

Now I want to export this project using qgis2web so that the map shows the pictures just like it does in QGIS.
What will I have to do??


texts: Launching multiple instances of QGIS (Mac)?

texts: Often I need to work on two qgis projects at the same time. How can I launch two separate instances of QGIS?


texts: Label Rotation in GeoServer SLD based on CRS

texts:
I'm styling a series of vector data in Polar Stereographic Projection (specifically EPSG 3031) using GeoServer and SLDs.  These maps can be rotated E/W, so if the labels are fixed they appear upside-down in places. (Example here: SCAR ADD (hold shift and drag to rotate)).

Is there a function that I could apply to the Rotation feature in SLDs to calculate the rotation of the labels based on the CRS rotation/meridian of the map?


texts: extract data from several images to different shapefile columns using GRASS

texts:
I am now working in grass and I want to extract raster values from several images to a shapefile. I have already created the columns in the .shp where the raster values will be stored but I am facing difficulties on how will I make the loop...

What I want is to extract the values of image_1 to column_a, values of image_2 to column_b etc etc... The problem that I have is that I can't make the loop for two variables (raster image and shapefile column) together... It must be easy but I am new in grass and bash in general as you must have already understood.

The command that I used is "v.what.rast"

Thank you!


texts: What is the need of "onExtentHistoryChange" in map navigation toolbar

texts:
In map navigation toolbar what is the need of "onExtentHistoryChange" in on(navToolbar, "onExtentHistoryChange", extentHistoryChangeHandler);

<!DOCTYPE html>

<meta http-equiv="Content-Type" content="text/html; charset=utf-8">

<meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
<title>Navigation toolbar</title>

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

Zoom In

Zoom Out

Full Extent

Prev Extent

Next Extent

Pan

Deactivate


texts: Create rectangular polygon with dimensions in km centered on points in point layer

texts: Given a point layer in QGIS 3.03, I am attempting to create a second layer of square polygons with a specified length of x kilometers whose centers are the points in the first layer.  What would be the best way of accomplishing this?


texts: Can I Turn Off or Remove Map Layers in QGIS Print Composer

texts: The title pretty much says it all.  I have a base map in QGIS that I am using with several Print Composers.  In one of the Print Composers, I don't want one of the layers to show up.  Is there any way to accomplish this?


texts: Which training courses in geostatistics can you recommend?

texts: I'm looking for excellent geostatistics courses all over the world, including online courses. In particular, I would like to learn more about the kriging method. Can anyone give me some suggestions?


texts: WFS vector layer not rendering with OpenLayers?

texts:
I'm having problems trying to load a Geoserver/WFS vector layer using OpenLayers. I can see GML data loading in Firebug and also I can see the layer name is in the layer switcher but the the data is not rendering. I've added my code below. I've seen similar questions on here that have helped me a lot with this so far but I can't see anything else to make it work.

Can anyone see anything I've missed or done wrong?

    wfs1 = new OpenLayers.Layer.Vector("CamEO Boundary", {
        strategies : [new OpenLayers.Strategy.Fixed()],
        protocol : new OpenLayers.Protocol.WFS({
        version : "1.0.0",
        url : "http://domain.co.uk/geoserver/wfs",
        featurePrefix : "riversdata",
        featureType : "cameo",
        featureNS : "http://domain.co.uk/riversdata",
        srsName: "EPSG:27700",
        geometryName : "geom"
        }),
        format : new OpenLayers.Format.GeoJSON({
        'internalProjection': new OpenLayers.Projection("EPSG:27700"),
        'externalProjection': new OpenLayers.Projection("EPSG:900913")
      }),
        styleMap: new OpenLayers.StyleMap({
        strokeWidth: 3,
        strokeColor: "#333333" })

      });

map.addLayer(wfs1);


texts:
I´m using QGIS Lyon 2.12 and am trying to make a map in the print composer with the scale 1:25000 in A4 size. I have drawn the map layer, but the map only covers a small part of the A4 space in the correct scale. I´d like to resize the layer in print composer so that it covers more space, while retaining scale and aspect ratio. Is there a way to do this?

If not, is it possible to crop the page while keeping A4 proportions? Basically, cut away all the white areas.


texts: What is stored in Options|Advanced|PostgreSQL|connections|dbname|keys group in QGIS 3?

texts:
Since QGIS 3 in the Options|Advanced|PostgreSQL path referred to in the title there is a new group 'keys' as can be seen in this screenshot:

In old QGIS 2.X:

What kind of information is stored in this 'keys' group and what do the Types / Values mean?


texts: Delete Turned Off (Unchecked) Layers from Arcmap with Python

texts:
I have about 100 figures that have been updated countless times and a lot of the old data was left in each .mxd but was only turned off (unchecked) and not removed.  I would like to remove these layers from each .mxd but keep the actual shapefile/features in an achieve folder.

I have found a similar script that removes layers by layer name here:
Delete layer from ArcMap using Python

but it does not do exactly what I want it to do.

Does anyone have a script already written that does what I am trying to do?


texts: How to import ESRI Ascii format into GRASS?

texts:
I need to take a set of datafiles, generated by our sensor system that reads GPR (Ground penetrating Radar) deviations over a surface as we drive along.  I have discrete readings from a set of sensors and can accurately place them relative to the GPS unit spitting out location values once a second, along with track and speed in meters/second.  I want to create an Esri ASCII data file for each second, with each row representing the sensor values along the rail, with one row representing all 16 sensors (spaced equally apart), read every 1/25th of a second.  I can figure out the offsets to accurately place the sensors, in meters, relative to the GPS unit and as each sensor is read.  It is relatively easy to then figure out the location for each reading along the track the GPS unit is measuring in Meters/second, thus arriving at the separation, in meters, between rows of generated sensor data.

What I am having a hard time doing is converting this information into ESRI ASCII format, specifically the CELL size.  My sensor rack is about 4 meters long and the GPS unit is centrally located with a small offset forward of the sensor array.  How can I determine the cell size for my many rows of sensor data?  I will locate this data with the GPS output, and have no problem calculating the boundaries for the one second time period.  What has me stumped is the cell size.

Example data:
Sensors 1 - 16
232.2 233. 244.4 256.8 260.1 ..... 211.5
next row 0.1 meter later
next row 0.1 meter later
25 rows in all

Sensors separated 0.25 meter from each other

Any ideas folks?  What does that pesky cell size have to be and more importantly, how did you figure out that number?


texts: Unsupported WFS request in MapServer at Insert and Update

texts:
I'm created a WFS with MapServer (ms4w) and SQL Server DB. The wfs:GetFeature works fine. But I cann't get the wfs:insert/wfs:update to work, they return:

    <?xml version="1.0" encoding="ISO-8859-1"?>
    <ows:ExceptionReport xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows" version="1.1.0" language="en-US" xsi:schemaLocation="http://www.opengis.net/ows http://schemas.opengis.net/ows/1.0.0/owsExceptionReport.xsd">
      <ows:Exception exceptionCode="InvalidRequest" locator="request">
        <ows:ExceptionText>msWFSParseRequest(): WFS server error. Unsupported WFS request.</ows:ExceptionText>
      </ows:Exception>
    </ows:ExceptionReport>

My table:

    CREATE TABLE Plots
    (
    Id INT PRIMARY KEY IDENTITY,
    Title VARCHAR(30),
    Plot GEOMETRY
    )

and created a map file:

    MAP
      NAME           "crop"
      STATUS         ON
      CONFIG "PROJ_LIB" "C:\ms4w\proj\nad"
      WEB
        IMAGEPATH "c:\ms4w\tmp\ms_tmp"
        IMAGEURL "/ms_tmp/"
        METADATA
          "wfs_title"            "Crop plots" ## REQUIRED
          "wfs_onlineresource"   "http://localhost/mapserver/cgi-bin/mapserv.exe?map=/ms4w/apps/crop/crop.map&" ## Recommended
          "wfs_abstract"         "Crop plots" ## Recommended
          "wfs_enable_request"   "*"  # necessary
          "wfs_namespace_uri"    "http://www.cropvision.com"
          "wfs_namespace_prefix" "crop"
        END
      END
      PROJECTION
        "init=epsg:28992"
      END

      LAYER
        NAME         "plots"
        STATUS       ON
        TYPE         POLYGON
        METADATA
          "wfs_title"           "Plots"
          "gml_include_items"   "all"
          "wfs_enable_request"  "*"
          "gml_geometries"      "plot"
          "gml_plot_type"       "multipolygon"
          "gml_featureid"       "Title"
        END

        PROJECTION
            "init=epsg:28992"
        END
        PROCESSING 'CLOSE_CONNECTION=DEFER'
        CONNECTIONTYPE PLUGIN
        PLUGIN "C:\ms4w\Apache\specialplugins\msplugin_mssql2008.dll"
        CONNECTION "server=myServer;database=myDB;uid=myUID;pwd=myPWD"
        DATA "Plot FROM dbo.Plots USING UNIQUE id"
     END
    END


texts: I have a GeoJson feature that is rendered with leaflet and proj4(ESPG:31467). Now I want to rotate that feature so that it NW corner points to NE, e.g. 90 deg rotation. Is there a way to do this? Or can you give me hints how to address this issue and develop a LL plugin for such functionality.


texts: I don't understand the SAGA-module "NORMALISATION". In SAGA Gis I just can define one MIN and one MAX value. In QGis I have to define 2 of both. How to handle that exactly?


texts: Writing comments in QGIS Graphical Modeler

texts:
I want to explain my model step by step so anyone can easely understand what's going on.
So, is it possible to write comments in the QGIS's Graphical Modeler?


texts: Unzipping OSM2VectorTiles after extracing with MBUtil

texts:
After extracting the tiles from an OSM2VectorTiles set, I have a proper folder structure with .pbf-files. However, OpenLayers cannot read those files.
I believe this is because they are zipped. So, whats the proper way to extract those into unzipped usable .pbf-files?


texts: Calculating length of field using Arcade?

texts:
I have data with this format: 4024-457 and I need to separate the left and right portions from the '-'. The data is not always in that proportion (4 numbers left and 3 to the right) so I need to calculate the position of '-' to extract both sides.

In the case of left numbers I used this expression:

Left($feature.R_ROL,Find('-', $feature.R_ROL))

But i don't find the way to extract the right portion of the main data. In ArcMap I used a VBscript expresion like this.

Right( [R_ROL] , Len( [R_ROL])- (InStr( [R_ROL] , "-")))

I can't find anything like the Len function to calculate the total lenght of the field to recreate something like the above. I'm still new to Arcade.


texts: Raster: how to locate cells having nodata

texts:
I am new to oracle spatial. I am using SDO_GEOR package's getNoData prcedure. the results are populated in following format. MDSYS.SDO_RANGE(MDSYS.SDO_RANGE(-32768,null))

My questions are:

Is result an array object / just a some form specification ?
if an object how can I use / consume it further say to generate a list/ array of cell coordinates.

My main intention is to locate the cells having nodata from the oracle( original raster file had nodata / voids cells).

I am not a GIS/ spatial pro.

Query:
select SDO_GEOR.getNODATA(georaster, 0 ) from <> where id=<>


texts: How to flip a polyline based on a separate polyline feature class being left or right to it

texts:
Can't go into detail about what these features actually are but I will just refer to them as the green line and the black line.  The black lines are not able to ever overlap the blue polygons.  Black lines have a specific orientation depending on which side the closest green line is.

What I'm trying to do is create a script which would be able to determine if the closest green line is either right or left to the black line, then determine if the black line is oriented correctly or not based on the rules I drew out within the attached image.  If it isn't then it would flip the black line.

I've used the FeatureVerticesToPoints_management tool and SplitLine_management tool to split both the polylines into individual parts and retrieve midpoint X and Y coordinates as well as individual vertex X and Y coordinates of both lines.

I've been using the formula (Bx - Ax) * (Cy - Ay) - (By - Ay) * (Cx - Ax) where A and B are line vertices of the green line and C is the midpoint vertex of the black line.  This has created a positive or negative value signifying which side point C is on compared to the line but I do not know how to setup rules for how the black line should be orientated afterwards. I believe the answer lies within checking the black line's azimuth.

My original idea was that using that formula I could determine which side the black line was on, then check if the azimuth of the black line is within 180deg of the the azimuth of the green line, and if it wasn't, it would flip the black line.  But since the green line can be orientated in it's opposite direction, and be acceptable, the 180deg rule turns into a 360deg rule and the black line can then go any direction it wants.

I've been calculating the azimuth of each line using:

def CalculaAzimuth(linea):
if (hasattr(linea,'type') and linea.type == 'polyline'):
    xf = linea.firstPoint.X
    yf = linea.firstPoint.Y
    xl = linea.lastPoint.X
    yl = linea.lastPoint.Y
    dX = xl - xf
    dY = yl - yf
    PI = math.pi
    Azimuth = 0 #Default case, dX = 0 and dY >= 0
    if dX > 0:
        Azimuth = 90 - math.atan( dY / dX ) * 180 / PI
    elif dX < 0:
        Azimuth = 270 - math.atan( dY / dX )* 180 / PI
    elif dY < 0:
        Azimuth = 180
    return Azimuth
else:
    return False


texts:
I can define a custom CRS in QGIS (v2.18.13, Windows 8) and conversion works fine and no error is reported when I hit OK.

proj4 is +proj=lcc +lat_1=49 +lat_2=77 +lat_0=0 +lon_0=-95 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

But when I go to use this CRS, there are no user-defined projections in the list.

A similar issue was raised here and here, but neither answer fixed my problem. May be a permissions issue, as the problem is not limited to the specific proj4 that I use, but I've enabled write access on all QGIS folder, still with no luck.

I'm experienced with GIS applications, but new to QGIS.


texts: How can I Use Python in ArcGIS to Manipulate the Underlying Data of Selected Polygons

texts:
I am relatively new to python in ArcGIS.

I have a shapefile loaded in ArcGis (v 10.0) and will select two adjacent polygons using the selection tool. Once I select my two polygons, I want to use Python to examine the underlying records for those selected polygons, pull out some data from certain fields in the underlying attribute table, do some simple math using those data, and then write that result back out to the same attribute table (in an existing field  - i.e. I will overwrite data in that field).

Here's the code I have so far:

import arcpy
#this is the name of the shapefile that's loaded in ArcMap
lyr_nm = "layer name"

#Get the current map layers
mxd = arcpy.mapping.MapDocument("CURRENT")
lyr = arcpy.mapping.ListLayers(mxd, lyr_nm)[0]
#this gets the layer but I need a feature class to examine the underlying records

#set field variable names
acresField = "ACRES"
tcovField = "TreeCover"
sCovField = "ShrubCover"
sourceField = "Source"

try:

    #Check for exactly 2 selected features.
    desc=arcpy.Describe(lyr)
    if len(desc.FIDSet) != 2:
        arcpy.AddMessage("Please select only two polygons to merge!")
        sys.exit()

    napaAc = 0
    knoxAc = 0
    totalAc = 0
    napaTree = 0
    napaShrub = 0
    knoxTree = 0
    knoxShrub = 0
    treePct = 0
    shrubPct = 0

    #somewhere here I need to define a feature class based on the layer above,
    #and the feature class must only have those two
    #selected records in its result, otherwise don't this will work

    rows = arcpy.SearchCursor(fc)

    for row in rows:
        if row.getValue(sourceField) == 0:
            napaAc = row.getValue(acresField)
            totalAc = totalAc + napaAc
            napaTree = row.getValue(tcovField)
            napaShrub = row.getValue(scovField)
        else:
            knoxAc = row.getValue(acresField)
            totalAc = totalAc + knoxAc
            knoxTree = row.getValue(tcovField)
            knoxShrub = row.getValue(scovField)

    #calculate values to write back out to Knoxville record
    shrubPct = ((napaAc/totalAc)*napaShrub) + ((knoxAc/totalAc)*knoxShrub)
    treePct = ((napaAc/totalAc)*napaTree) + ((knoxAc/totalAc)*knoxTree)

    rows = arcpy.UpdateCursor(fc)

    for row in rows:
        if row.getValue(sourceField)==1:
            row.setValue(tCovField,treePct)
            row.setValue(sCovField,shrubPct)
            rows.UpdateRow(row)

except:

     # Report an error messages
     arcpy.AddError("Could not complete the polygon join!")

     # Report any error messages
     arcpy.AddMessage(arcpy.GetMessages())


texts: Add-in; Failure to get code to execute when toolbar button pressed

texts:
I had written several VBA routines for ArcGIS up through v9.3 and have since upgraded to 10. Since VBA will be a thing of the past in the next version, I downloaded and installed VB Express 2008 and successfully (I think) migrated the VBA to VB.Net 2008. I saved out each of the routines to their respective *.vb files and ran the Add-In wizard for a new project, adding each routine as a button added to a toolbar.

Everything compiles and I get my toolbar to show up in ArcMap, but when I press any of the buttons, they grey out and nothing executes. Similarly, the breakpoints I have inserted go hollow with a yellow triangle and when I hover over the breakpoint, it states that "The breakpoint will not currently be hit. No symbols have been loaded for this document."

I think that somehow I am missing the magic of how to tie the XML stuff (which I understand controls the toolbar and buttons) to the VB.Net code (which perform the actions). Is there any VB code that links the code back to the XML so it functions?

Any insight is greatly appreciated.

Thanks,

Dale Bridgford

P.S (added later) To provide a bit more background, I am enclosing a portion of the code.  This is typical of what I have cludged together from various posts.  All variations that I have tried are equally inoperative.

Again, any insight as to where I have gone astray would be greatly appreciated.

Thanks,
Dale Bridgford

`
Imports ESRI.ArcGIS.Geometry

Namespace ERMDTools

Public Class btnAddHyperlink
    Inherits ESRI.ArcGIS.Desktop.AddIns.Button
    Public Sub New()

    End Sub

    Protected Overrides Sub OnClick()
        Dim pFeature As IFeature, pLayer As ILayer
        Dim pMap As IMap
        Dim pSRE As SpatialReferenceEnvironment
        Dim pQF As IQueryFilter
        Dim pMSR As ISpatialReference
        Dim pGeoDataset As IGeoDataset
        Dim pMxDoc As ESRI.ArcGIS.ArcMapUI.IMxDocument
        pMxDoc = Nothing
        pMxDoc = m_App.Document
        pMap = pMxDoc.FocusMap
        pMSR = pMap.SpatialReference `


texts: Raster attribute table not showing up after exporting from R

texts:
I exported a categorical raster from R and loaded it into ArcMap to visualize it. However, the attribute table is greyed out.

Code In R:
I first tell R that the values are factors by ratifying it.

lulc_1985.resam <- ratify(lulc_1985.resam)
rat1985 <- data.frame(
  ID= 1:11,
  LandCover = c("Deciduous","Cropland", "Settlement","Mixed Forest",
                "Scrubland","Barren Land","Wetlands","Plantations",
                "Coastal","Grassland","Evergreen")
)
levels(lulc_1985.resam)<- rat1985

The output of the above:

> lulc_1985.resam
class       : RasterLayer
dimensions  : 203, 192, 38976  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 628560.6, 820560.6, 1084984, 1287984  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : layer
values      : 1, 11  (min, max)
attributes  :
       ID LandCover
 from:  1 Deciduous
 to  : 11 Evergreen

As you can see above, the attributes exist here. Now I export the same out of R:

writeRaster(lulc_1985.resam,"E:\\lulc_1985.tif", format="GTiff")

The raster is now visualized in ArcMap, where the attribute table is greyed out.

How do I retrieve this attribute table?

The exported .tif is just 13kb in size.


texts: Remove duplicate polygons in a single layer using ogr2ogr?

texts: I have a spatialite database file containing buildings within a country as multipolygons. While trying to use QGIS to delete these duplicate features, QGIS keeps on crashing due to the millions of buildings within the layer. So I am trying to find an command line solution using ogr2ogr but cannot find any info. What would be the syntax to do such a task?


texts:
I'm making a web application with leaflet and OSM. It's all going well, but now I need a little help.

I'm searching for a good routing engine, that's capable of creating the most probably taken route based on "checkarea":

I've many points logged every N minutes, and they are not really precise, but their offset should be 400m, average of 200m. I need to know what route has been taken.

I'll explain better:

I've this start A and finish B point. The black line it's like a street map, blue dots are points I calculated and blue circle are areas where it should route. The yellow path it's what the engine should calculate.

I've starting modifying Pyroutelib2, but maybe I'm just reinventing the wheel.


texts: Calculating spaces between roads from PostGIS database?

texts:
Is it possible to calculate the spaces between the roads using a select statement in a PostGIS/PostgreSQL Database?

I tried using some combinations of the predefined functions as ST_Intersects & ST_Difference but did not succeed yet.


texts: Calculate area covered by a particular raster layer within polygons in QGIS

texts:
I am trying to extract information on vegetation cover within a set radius of point locations (my sites). I have succeeded in creating a new vector layer of circles with a 5 km radius around my sites.

What I want to do is calculate the amount of vegetation cover within each of these circles (and therefore, within a 5 km radius of my sites), which I have as a raster layer. I would like to extract this either as a direct measurement of land cover (e.g. metres squared), or even just a percentage of the 5 km radius circle/polygon that is covered by the raster layer of vegetation (because I could then calculate the land cover from that).

I am aware that the LeCos plugin can somehow do this, but I have not succeeded in installing this (subject of another question I have posted) and was hoping there might be another way. I would be open to trying this in R if someone can point me in the right direction ( I am a user of R for statistical analysis but haven't yet tried it for spatial analysis).


texts: Always syntax error using query builder - QGIS3

texts:
I have a problem with QGIS 3.4.12. As you can see in the image below, it's not possible to execute queries using query builder, because there's always a syntax error. This only happens using this layer ("benigeo"). How can I solve this issue?

EDIT: I need to select the list of some attributes, for this reason I don't need a WHERE clause. To be more clear, I attach another screenshot where you can see what I need. In this layer ("waterways") I'm able to SELECT the waterways type list. Instead, in the other layer I have always that error.


texts: Load data from URL to update Billboard position Cesium

texts:
I am trying to update the position of a billboard every second in Cesium. The issue is I cannot pass the updating points to the billboard coordinates. I am using jQuery to load the data. This is what I have

//getting values every second
window.setInterval(function nameSat()
        {
         $.getJSON('http://10.175.129.10/meo_coordinates_heo.json', function(data)
            {
                 lat1 =`${data[0].Latitude  }`
                 lon1 =`${data[0].Longitude }`
                 alt1 =`${data[0].Altitude *1000 }`
                 active1 =`${data[0].Active }`

                var lat = parseInt(lat1);
                var lon = parseInt(lon1);
                var alt = parseInt(alt1);
            })

    }, 1000);

//i want to put the updating coordinates here
function addBillboard() {

    Sandcastle.declare(addBillboard);
viewer.entities.add({
        position : Cesium.Cartesian3.fromDegrees(lon180, 40,8000000),
        billboard :{
         scale : 0.03,
            image : '../images/sat.png'
        }
    });
}


texts: What is the correct syntax for the geometry object's "intersect" method in ArcPy?

texts:
I'm building a script that utilizes the intersect method of the geometry class so that I can run an intersect analysis on individual polygons within a feature class. As stated here,
 "intersect (other, dimension)" requires two parameters, but there is no explanation as to what kind of string they are expecting for "dimension". What should I pass to the method for that parameter? I can't find any documentation to tell me :(


texts: Disconnect Polygons connected to each other by roads?

texts:
I am working with a land cover polygon.  I have a "land" class and a "road" class.  I am trying to disconnect islands in this polygon from each other.  They are connected by roads. Can anyone think of a way to disconnect the islands, without removing the roads from the polygon completely. I do not want to remove them completely, as this would create holes in the islands where a road is. I am trying to avoid having to manually cut the roads out in Editor.  Below is a screen shot of the situation I am talking about. The data set is for a large area. and a couple thousand islands.  I am working in Arcmap 10.4.1

The tan areas are islands (land) the black is the roads connecting the islands.


texts: Domains tab missing in geodatabase properties?

texts:
I'm fairly new to ArcGIS Desktop. I started working on a feature layer in a new geodatabase yesterday, adding attribute fields and creating coded domains for them.

However today when I opened it up again, the Domains tab is missing. ie. it doesn't show up in Catalog --> mygdb.gdb -> [right click] -> Properties. Yesterday there were only 2 tabs in Properties - General & Domains. Today there's General, Sharing, Security, Previous Versions, Customize. But no Domains.

Any idea why this is or how to get it back?

I know I can go into Domains from the Data Management toolbox to create them but it's quicker and much more convenient to do it from the gdb's Properties.


texts: Query after parametric geometries in a GDB

texts:
When you create an circle or ellipse in a GDB it is rendered as a parametric geometry storing x, y, radius, etc. and not as a polygon with nodes.

http://desktop.arcgis.com/en/arcmap/10.3/manage-data/using-sql-with-gdbs/parametric-circles-and-ellipses.htm

I have a layer with both polygons and parametric geometries (circles/ellipses). How do I make a query in ArcGIS 10.3 that returns only the parametric geometries and not the 'normal' node based polygons?


texts: Generating UUID(unique ID) at joining points

texts:
I have road network (lines) and points (Address) . With the help of closestSegmentWithContext function, I join points to nearest lines in QGIS. Is there a way to to generate a unique ID at every joining point as the joining points are drawn manually in below image.
What I have researched so far:

Found a possible solution  here, but it didnt worked for joining points.
I come accross a useful plugin "AUTOFIELD" in anwer here, but this plugin is no more avaible in QGIS 3.XX

The code for joining line segements to nearest points is given below and shape files are available here:
roadlyr = QgsProject.instance().mapLayersByName('TR_ROAD')[0]
addrlyr = QgsProject.instance().mapLayersByName('ADDRESS')[0]

#List all features
adresses = [f for f in addrlyr.getFeatures()]
roads = [f for f in roadlyr.getFeatures()]

#Create empty vector layer
vl = QgsVectorLayer("LineString?crs={}&index=yes".format(roadlyr.crs().authid()), "myLayer", "memory")
provider = vl.dataProvider()

for adress in adresses:
    a = min([road.geometry().closestSegmentWithContext(adress.geometry().asPoint()) for road in roads], key=lambda x: x[0]) #For each road find closest address
    gLine = QgsGeometry.fromPolyline([QgsPoint(adress.geometry().asPoint()), QgsPoint(a[1])])
    f = QgsFeature()
    f.setGeometry(gLine)
    provider.addFeature(f)

QgsProject.instance().addMapLayer(vl)


texts: Creating DEM from multiple .las files using LAStools in ArcMap?

texts:
I am attempting to create a DEM for a raster large dataset, with over 500 .las files. Within the las2dem tool, within the LAStools toolbox there is not an option to select an entire folder, nor can you do so manually.

Thus, I tried creating a DEM for each file; however, an output was not generated. I saved both the file location and directory, but got no output. How do I use this tool?

Also I attempted to merge some of my .las files together using the lasmerge tool and got no output. I also ensured the file location an directory was saved properly.