texts: Extracting NDVI scores over time for defined field boundaries

texts:
I am trying to do something similar to this post with the dataset "users/ryanhuffman/Training_Dataset_Trunc" and keep getting hung up on an error with the dictionary section of the code. For simplicity I renamed my column headers and kept variable names the same to streamline retro-fitting code.
Link below for context.
var Akdarya = ee.FeatureCollection("users/ryanhuffman/Training_Dataset_Trunc")

Map.centerObject(Akdarya)
Map.addLayer(Akdarya)

var S2 = ee.ImageCollection('COPERNICUS/S2')
  .filterDate('2018-05-01', '2018-09-30')
  .filterBounds(Akdarya)
  .filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 20)

var maskcloud = function(image) {
  var QA60 = image.select(['QA60']);
  return image.updateMask(QA60.lt(1))
}

var NDVI = function(image) {
  return image.expression(
    '(NIR - RED) / (NIR + RED)', {
      'NIR': image.select('B8'),
      'RED': image.select('B4'),
    }).rename('NDVI').copyProperties(image, image.propertyNames());
}

var NDVIcol = S2.map(maskcloud)
  .map(NDVI)

var ts = NDVIcol.map(calculateFarmersNDVI)

function calculateFarmersNDVI(image) {
  var properties = ee.Dictionary(
    Akdarya.iterate(function(region, acc) {
      return ee.Dictionary(acc)
        .combine(
          ndviForArea(image, ee.Feature(region))
        )
    }, ee.Dictionary())
  )
  return ee.Feature(null, properties
    .set('Year', image.date().get('year'))
    .set('Month', image.date().get('month'))
    .set('Day', image.date().get('day'))
  )
}

function ndviForArea(image, region) {
  var ndvi = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region.geometry(),
    scale: 30
  }).get('NDVI')

  var farmers = region.getString('Farmers')

  var ndviOrBlank = ee.Algorithms.If(ndvi, ndvi, '')
  return ee.Dictionary().set(farmers, ndviOrBlank)
}

var listOfFarmers = ee.List(Akdarya.aggregate_array('Farmers'))
  .map(function (farmers) {
    return ee.String(farmers)
  })
ee.List(['Year', 'Month', 'Day'])
  .cat(listOfFarmers)

  .evaluate(function (selectors) {
    Export.table.toDrive({
      collection: ts,
      description: 'NDVI_test3',
      fileFormat: 'CSV',
      selectors: selectors
    })
  })

Exporting NDVI result data as a table using Earth Engine


texts: How to overlay point co-ordinates on tmap shapefile?

texts:
So far, I've managed to visually present data for different socio-economic factors at a district level in Pakistan with the help of the tmap package in RStudio to create the following output.

At the same time, I have used the ggmap package to plot the geographic co-ordinates of various towers for different Mobile Network Operators in Pakistan (roughly 13000 towers). The output is as follows.

My next requirement is to overlay the second plot over the first, essentially showing the towers on top of the district wise electricity scores.

Any ideas on which packages I could use (as well as sample code) in order to achieve this?


texts: Using Select Layer by Location in ModelBuilder?

texts:
I'm trying to export a selection ( "select by location") using Copy features.

The selection is from the "Layer1".

It apparently executes the selection, but the copy feature is made by the original layer and not by the selection.

I am using NEW_SELECTION in "selection type" of the "Select Layer by location" menu.

In the "Copy Features" the input feature is equal to the original file name.

Any suggestions?


texts:
I have a tif with the following projection:

'+proj=stere +lat_0=1 +lon_0=2 +k=1 +x_0=6 +y_0=7 +datum=WGS84 +units=m +no_defs'

Geoserver doesn't accept these tifs and QGIS is unable to reproject them.

How do I know if this projection is ok and what do I do to correct it if not?

I add the gdalinfo output for reference:

Driver: GTiff/GeoTIFF
Files: xxxxx.tif
Size is 2751, 826
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Stereographic"],
    PARAMETER["latitude_of_origin",1],
    PARAMETER["central_meridian",2],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",6],
    PARAMETER["false_northing",7],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-18.500000000000000,46.500000000000000)
Pixel Size = (0.019992729916394,-0.019975786924939)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -18.5000000,  46.5000000) (  1d59'59.21"E,  1d 0' 1.29"N)
Lower Left  ( -18.5000000,  30.0000000) (  1d59'59.21"E,  1d 0' 0.75"N)
Upper Right (  36.5000000,  46.5000000) (  2d 0' 0.99"E,  1d 0' 1.29"N)
Lower Right (  36.5000000,  30.0000000) (  2d 0' 0.99"E,  1d 0' 0.75"N)
Center      (   9.0000000,  38.2500000) (  2d 0' 0.10"E,  1d 0' 1.02"N)
Band 1 Block=2751x1 Type=Byte, ColorInterp=Red
Band 2 Block=2751x1 Type=Byte, ColorInterp=Green
Band 3 Block=2751x1 Type=Byte, ColorInterp=Blue


texts: Inserting lake in height map using QGIS?

texts:
I am trying to merge two height maps in QGIS.

I have DEMs of the lake of constance in germany. It is a result of a goverment project and the data is available in a asc format. I want to 3dprint a landscape relief model. Also i need to split that model into pieces (200x200mm) to print it in blocks for a total length of 1000cm.

There is no problem to import the data in Qgis, also it is no problem to split the height maps into pieces for 3d printing. The fun begins when i try to merge the two DEMs into one. The DEM of the surrounding area ends where the lake begins.

Is it possible to "insert" the lake into the image of the surrounding area and recolor it?

The deepest point in the new DEM should be the darkest in the whole new DEM. At the end i need a single height map. --> Target: DEM of the whole area wir a dried of lake.

Shows the lake in the middle plus the surrounding area.

Shows only the surounding area


texts: What Z-factor to set in Curvature (spatial analyst)?

texts:
I'm making a slope curvature map, my DTM is 5*5 meters cells, in a meter grid projection.
The elevation is in meters asl.
When using the Curvature tool for different units of measure, the Z-factor has to be set.

  If the x,y units and z units are in different units of measure, the
  z-factor must be set to the appropriate factor, or the results will be
  incorrect. For example, if your z units are feet and your x,y units
  are meters, you would use a z-factor of 0.3048 to convert your z units
  from feet to meters (1 foot = 0.3048 meter).

Yet I'm not exactly sure if I need to set a Z-factor of 5, or use the default of 1, because in the "how curvature works" page, the illustration shows that the value for "L", is one raster cell (in my case, 5 meters)

What is your opinion?


texts:
When I select a geometry on map, the method map.infoWindow.getSelectedFeature() always returns "undefined". It looks like it is not working.
This behavior is with InfoWindow and also with Popup.
Also not working event map.infoWindow.on("selection-change",....

Code:

map = new Map("mapDiv", {});
map.infoWindow.on("selection-change", function ()
{
   var graphic = map.infoWindow.getSelectedFeature();
   // graphics is always "undefined"
}

How to read currently selected feature?


texts: What does ArcGIS create in a File Geodatabase (.gdb) when rasters are added to a mosaic dataset?

texts:
Our team frequently creates mosaic datasets, and adds rasters to them. We keep older datasets (and rasters) for archival purposes.

Now that we're running low on disk storage, we'd like to remove some of the older rasters. But we're not sure whether doing so breaks existing services. To that end, we'd like to know what ArcGIS creates when rasters are added to a mosaic dataset. Does ArcGIS replicate the rasters into the GDB folder? Or does it simply create metadata that points to the rasters at their source?


texts: Disable feature popup when creating new simple marker

texts:
I have a leaflet map with some loaded polygon layers.
When using the Simple Markers plugin, i would like to create a new marker on top of an already existing polygon feature, but i can only click the feature (which brings up the popup with information) and not place the marker on top of the feature. It works fine on an empty piece of the map.

Is it possible to disable the feature popup when the Simple Marker function is running?

L.Control.SimpleMarkers = L.Control.extend({
    options: {
        position: 'topleft'
    },

    onAdd: function () {
        var marker_container = L.DomUtil.create('div', 'marker_controls');
        var add_marker_div = L.DomUtil.create('div', 'add_marker_control', marker_container);
        add_marker_div.id = 'add_marker_div_id';
        add_marker_div.title = 'Add a new Note';
        L.DomEvent.addListener(add_marker_div, 'click', L.DomEvent.stopPropagation)
            .addListener(add_marker_div, 'click', L.DomEvent.preventDefault)
            .addListener(add_marker_div, 'click', (function () { this.enterAddMarkerMode() }).bind(this));

        return marker_container;
    },

    enterAddMarkerMode: function () {
        if (markerList !== '') {
            for (var marker = 0; marker < markerList.length; marker++) {
                if (typeof(markerList[marker]) !== 'undefined') {
                    markerList[marker].removeEventListener('click', this.onMarkerClickDelete);
                }
            }
        }
        document.getElementById('map').style.cursor = 'crosshair';
        document.getElementById('add_marker_div_id').style.backgroundColor = '#b0c4de';
        map.addEventListener('click', this.onMapClickAddMarker);
    },

    onMapClickAddMarker: function (e) {
        map.removeEventListener('click');
        document.getElementById('map').style.cursor = 'auto';
        document.getElementById('add_marker_div_id').style.backgroundColor = '#FFF';

        var popupContent =  "You added a note on the map at " + e.latlng.toString();
        var the_popup = L.popup({maxWidth: 160, closeButton: false});
        the_popup.setContent(popupContent);

        var marker = L.marker(e.latlng);
        marker.addTo(map);
        marker.bindPopup(the_popup).openPopup();
        markerList.push(marker);

        return false;
    },
});
var markerList = [];

Edit:

I have added the features from a PHP file through the following code as a GeoJSON layer:

var exp_PolygonJSON = new L.geoJson(exp_Polygon,{
    onEachFeature: pop_Polygon,
    style: doStylePolygon
});

The popups are coming from:

function pop_Polygon(feature, layer) {
    var popupContent = '____Popup html code____';
    layer.bindPopup(popupContent);
}


texts: Clipping with gdal translate does not create a file (Python)

texts:
I have tried to clip a geotiff file using the code below:

from osgeo import gdal

miny=-10.616665787
maxy=38.558338147
minx=34.562256098
maxx=71.187259028

ds = gdal.Open(fp)
ds = gdal.Translate(out_tif, ds, projWin =  [maxx, miny, minx, maxy])
ds = None

For some reason no file is generated though, when I tried using different coordinates it worked perfectly:

from osgeo import gdal

ds = gdal.Open('original.tif')
ds = gdal.Translate('new.tif', ds, projWin = [-75.3, 5.5, -73.5, 3.7])
ds = None

Does anyone have an idea what the problem could be? This is the area I am working with:


texts: I am interested in logging NMEA data on a Trimble Nomad 900. Unfortunately the SatViewer records a maximum of 64,000 characters, too little for my application. I don't want to have to purchase Terrasync/Pathfinder because I already have tools to parse NMEA and those pieces of software are overkill for my intended use. Thanks in advance.


texts: QGIS - add Layer From URL (GetCapabilities)

texts:
I'm trying to add a layer to my QGIS project from this URL: http://quegis.labcd.unipi.it/cgi-bin/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities

Using a WMS connection, I'm able to open it in jpeg or PNG format, but I don't know how to retrieve the alphanumeric part, i.e. feature attributes, which is the interesting part for me. Am I making a mistake or is it impossible via WMS?


texts: Creating shortcut (button) for "Make this the only selectable layer" in Python 2.6?

texts:
I'm trying to create a script in Python 2.6 to make a button or a keyboard shortcut for "Make This The Only Selectable Layer". I came across the following VBA scripts on the net, but the problem is it's not working in ArcGIS 10. Can anyone help me with a Python script to do that?

Public Sub MakeLayerSelectable_Click()

  Dim pMxDocument As IMxDocument
  Dim pMap As IMap
  Dim pEnumLayer As IEnumLayer
  Dim pLayer As ILayer
  Dim pId As New UID
  Dim pFLayer As IFeatureLayer

  Set pMxDocument = Application.Document
  Set pMap = pMxDocument.FocusMap
  pId = "{6CA416B1-E160-11D2-9F4E-00C04F6BC78E}"
  Set pEnumLayer = pMap.Layers(pId, True)
  pEnumLayer.Reset
  Set pLayer = pEnumLayer.Next
  Do While Not pLayer Is Nothing
    If TypeOf pLayer Is IFeatureLayer Then
        Set pFLayer = pLayer
        pFLayer.Selectable = False
    End If
    Set pLayer = pEnumLayer.Next
  Loop

Set pLayer = pMxDocument.SelectedLayer

If TypeOf pLayer Is IFeatureLayer Then
  Set pFLayer = pLayer
    pFLayer.Selectable = True
End If

End Sub


texts:
I have a layer saved to spatialite which I want to query and save part of to csv.

Even the query result takes a very long time to load to the browser (2.5M rows) even if I load it without geometry.

Is there a way to use the SQL engine DB Manager has to save the query result to file?

I already tried the sqlite way with the example I found here.
and the COPY TO way PostgreSQL uses specified here.

Is there another way I am missing? I thought QGIS DB Manager uses an SQLite engine.

Using QGIS 3.2.1 on Windows.


texts: How can I edit name of lines of layer in ArcMap from Visual Studio C#? E.g., change STR_NAME = "Pushkin street" to STREET = "Pushkin street"?


texts: ArcGIS 10 - mosaic: creating 4 files (rasters)

texts:
using the Mosaic to new raster tool in order to merge two satellite images, ArcGIS 10 created 4 different files: the coloured image and image1, image2 and image3 (all BW), which are, I guess, three different bands.

How to merge them all to have only one file (raster)?
(for the purpose of storing the file in the database for the other users, to make it very simple)


texts:
How do I clip a shapefile with a raster file with gdalwarp?

Input: a shapefile
Cutter: a raster file
Output: a cut shapefile

I can only find the answer of cutting the raster file by a shapefile but not vice versa.


texts: Output extent environment is not respected by Select Layer By Attribute (ArcGIS 10.2)

texts:
I cannot get the arcpy.env.extent working for Select Layer By Attribute GP tool. I can verify that I get various values for coordinates when I have panned to various areas on the map. Copy Features GP tool does export only those feature that are located within the current visible map extent. According to the Select Layer By Attribute GP tool documentation:

  If an Extent environment is specified, or if a definition query is
  present on the Input Layer or Table View, only the features or rows
  within the extent or matching the definition query may be selected.

Has anyone succeeded selecting features that are within the extent only? When I run the code below, regardless of the map extent, all features that fulfill the query requirements are selected.

This is the code I execute in the Python window in ArcMap.

fc = "Water distribution network fittings" #layer name in the TOC
field = "ObjectId"
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]
print df.extent.XMin, df.extent.YMin, df.extent.XMax,df.extent.YMax

arcpy.env.extent = df.extent #also tried arcpy.Extent(XMin, YMin, XMax, YMax)
arcpy.CopyFeatures_management(fc, "in_memory\outputFeatures") #adds new layer correctly
arcpy.SelectLayerByAttribute_management(in_layer_or_view=fc,selection_type="NEW_SELECTION",where_clause = """ObjectId > 160""")


texts: Find intersection of multiple DSM's using Python and GDAL

texts:
I am quite a beginner in respect to GIS and hopefully someone can help me out with the problem I have.

I have multiple digital surface models [DSMs] of the same area with all necessary geo-information (6 in total). However they are not exactly the same as they have been acquired on different days (row and cols of the raster are not identical). I already managed to align the DSMs, however now I need to find the intersecting area that is part of all six DSMs because I want to do some processing with this (which I indent to do in python).
I found the following post:

http://sciience.tumblr.com/post/101722591382/finding-the-georeferenced-intersection-between-two

This methods works fine for two DSMs, however the extension to more areas is not straightforward (at least not for me at the moment).

I guess this is a rather standard beginner problem, however I only ever found answers regarding exactly two areas.


texts: Batch Mode Routing for Fastest Route Distance

texts:
I would like to run a batch process to derive the total driving distance between points, for a largish set of point-pairs.

e.g. I would have data something like:

Start_ID    Start_Lat    Start_Long    Dest_ID    Dest_Lat    Dest_Long
00000001    42.12345     -76.54321     10000000   41.12345    -75.54321

What I would like to get returned from this is the driving distance.

I have tried this using a few web services, but the latency is such that it would take too long for my purposes.  I have potentially millions of routes to compute, so a latency of a couple seconds quickly adds up to way too much! :)

I have seen some references to setting up an OSRM instance as well as other references to PostGIS with pgRouting.  I don't know these tools, so I am not sure if they would meet my requirements or if there is another tool that would be better.

So, what good solutions exist for generating total driving distance for many pairs of points, to execute in batch mode?


texts: Changing the default cache directory through ArcGIS Server server manager

texts:
I am attempting to change the default cache directory for ArcMap so that when I publish a service I can cache its tiles on an alternative drive to Local Disk due to insufficient memory issues. However when I tried to do this using ArcGIS Server Manager I encountered the following error:

  Failed to register the server directory 'MIPCache'. Cannot write to directory path 'Z:'. Please check that the location is valid and that the ArcGIS Server account has permissions to the location.

Please note that permission should be in place, that I can manually access the folder location and that the user carrying out this process is the same one as the user that installed ArcMap on the server.

I am at a roadblock and have no idea what may be wrong...


texts: QGIS/PostGIS Difference Geoprocessing with tolerance

texts:
This is probably an extended question from my another question (Increase spatial query performance (st_difference))

This is the background is this. I have a road link and a tunnel link.
From the road link, I need to erase the tunnel link out.
With builtin GIS geoprocessing, I can use these tools Erase tool in ArcMap, Difference toon in QGIS.

I used ArcMap erase tool, and it seemed to work fine.
with QGIS tool, some areas were fine but another areas did not work well.

I assume that major issue was the topology of two links were not the same. (It does look like but it was not).

I wonder why ArcGIS returned the answer that I want but QGIS/PostGIS didn't.
Also, I would like to have exact same answer of ArcGIS erase tool in QGIS (though I know it is difference but I wonder if there are any tricks...)

Here is one of the screenshots of ugly areas.
considering that the scale is negative 9832 (1: -9832) the topology error? is not viewable with ArcMap (I suppose)


texts: Why does pyproj give a different point location compared to ordnance survey when I convert from WGS84 to EPSG:27700 (BNG)

texts:
I get different results using 2 methods to convert co-ordinates from WGS84, epsg:4326 to BNG (British National Grid) epsg:27700. The differences in the plotted eastings and northings values are very small resulting in a difference of approximately 0.5m.
The main reason I am concerned is that when I plot the converted points in QGIS the results from using pyproj are in the same location as the non transformed points from WGS84. Converting the WGS84 points on the fly in QGIS plots the points in the same location as those I converted using the OS software.
I may need to use pyproj in the future and would like to ensure I have used the software correctly.
I am using pyproj and Grid InQuest II software from the Ordnance survey to run the conversions.
WGS84 Co-ordinate set

longitude = -3.209168460809744
Latitude = 55.950621342577172
Altitude = 116.378547668457031

pyproj result and code
x = 324588.97991822625 y = 673725.6892528223
trans1 = Transformer.from_crs("epsg:4326", "epsg:27700")
x1, y1, z1 = trans1.transform (lat, lon, alt)
print ("Transformer lon =   " +str(lon), " lat = " + str(lat) + " x = " + str(x1) + " y = " + str(y1))

Ordnance Survey result
Easting = 324589.0436663538 Northing = 673726.1910941075
OS conversion validation
I have 3 conversion methods from the Ordnance Survey Website

Grid InQuestII source code, used to create result above. https://www.ordnancesurvey.co.uk/business-government/tools-support/os-net/for-developers
Grid InQuestII software, a stand alone application, same co-ordinates but rounded to 2 decimal places
Online conversion page, same result to 3 decimal places. https://www.ordnancesurvey.co.uk/gps/transformation/

The OS website and code converts from from ETRS89, epsg:4937, not WGS84, epsg:4326. I did try using epsg:4937 in pyproj and got the same result as using epsg:4326.

OS comment on using ETRS89, EPSG:4937
In Europe, ETRS89 is a precise version of the better known WGS84 (EPSG
code 4326) coordinate reference system optimised for use in Europe;
however, for most purposes it can be considered equivalent to WGS84.
Specifically, the motion of the European continental plate is not
apparent in ETRS89, which allows a fixed relationship to be
established between this system and Ordnance Survey mapping coordinate
systems.

Rendering in QGIS
I plotted a .geojson file with my original WGS84 lat, long values to see how that location compared with the converted values. I load 2 .csv files, one with the OS conversion and one with the pyproj conversion.
Method 1
Project Properties - Datum Transformations EPSG:4326 to EPSG:27700 selected
WGS84 points and OS points in same locations.
pyproj points about 0.5m apart
Method 2
Project Properties - Datum Transformations blank
WGS84 and pyproj points in same location
This was the confusing result, I would have accepted a result where the OS conversion was slightly different from the pyproj conversion and the plotting of the WGS84 result was in a slightly different location.
Added images of project reference system and project properties


texts: Why QgsRasterCalculator result is incorrect in pyqgis?

texts:
I want to sum 4 raster layers by QgsRasterCalculator. Below code creates a raster layer as result, but pixel values in raster result is wrong. what is the error in this code?

from qgis.core import *
from PyQt4.QtCore import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

entries = []
rst1 = "home/Data/lyr1.tif"
rst2 = "home/Data/lyr2.tif"
rst3 = "home/Data/lyr3.tif"
rst4 = "home/Data/lyr4.tif"
myList = [rst1, rst2, rst3, rst4]

try:
    for i in range(0, len(myList)):
         raster = myList[i]
         readRst = processing.getObject(raster)
         ras1 = QgsRasterCalculatorEntry()
         ras1.raster = readRst
         ras1.ref = "ext_lyr@" + str(i+1)
         ras1.bandNumber = 1
         entries.append(ras1)

    formula = '(' + entries[0].ref + ' + ' + entries[1].ref + ' + ' + entries[2].ref + ' + ' + entries[3].ref + ')'
    readRst = QgsRasterLayer(myList[0])
    output = "home/Data/result.tif"
    calc = QgsRasterCalculator(formula, output, 'GTiff', readRst.extent(), readRst.width(), readRst.height(), entries)
    calc.processCalculation()

except:
    self.userWarning("Error in Create map", "Can not create Create map, Exit")


texts: Joining tables based on distance using SQL Server

texts:
I need to show the three closest hospitals to every house. I have a table of properties and a table of hospitals. Both are in MSSQL2012, have a geography column that is indexed and are in the same 4326 SRID.

I cant find any examples on how to perform ST_Distance across different tables, returning the top 3.

As per Evans suggestion i have tried the following but the query is not valid

SELECT *
FROM (
        SELECT * , row_number OVER (
        PARTITION BY [dbo].[houses].[geogtab]
        ORDER BY [dbo].[houses].[geogtab].STDistance([dbo].[hospitals].[geogtab]) ASC
      ) AS i
  FROM [dbo].[houses]
  JOIN [dbo].[hospitals]
       ON [dbo].[houses].[geogtab].STBuffer(100000).STIntersects([dbo].[hospitals].[geogtab])
       ORDER BY [dbo].[houses].[geogtab].STDistance([dbo].[hospitals].[geogtab]  ) ASC
) AS t
WHERE t.i <= 3;

I am unsure if this is formatting related and am unsure from reading this if the distance will be an output in the result table.


texts: Highlight polygon on hover and click in TileMill

texts: I have created a MapBox.js map with several layers of polygons each with hover/click interactivity to show attributes about each polygon.  I want the polygons to become highlighted upon a user hover/click with a darker border to help the user identify a polygon.  Something similar to this example.  What code would I need to add to make this happen? I'm assuming I need to create a CSS file to define the color and size of the outline but I don't know how to make it work with my JavaScript code.  My current source code is here.


texts: How to create a graphic point by xy value?

texts:
I have implemented this code to create point event layer by xy text file, however I'm wondering how to create a point graphic by xy value?

This script draws point graphic by mouse down.  When I steer away from the mouse down event and input my xy values manually it does not draw.
Using the Point class looks like it might also be a solution, however I'm not sure how to implement this code without getting a null reference error on point.setX line.

Dim point As New Point()
point.setX(-100.0)
point.setY(40)
Dim markerElement As New MarkerElement()
markerElement.setGeometry(point)

Dim map As IMap = mxdoc.FocusMap

Map.addElement(markerElement, 0)

Any suggestions on how to create a graphic point by xy value?

Thank You


texts:
Can somebody provide a good explanation/comparison between Info Window vs. Popup Window.
It seems both can be used with Feature layers. What are the pros and cons of either and when is more appropriate to use one than the other?


texts: Zoom into Attribute Table Itself in QGIS?

texts:
How do I zoom into the attribute table in QGIS? I need to see the cells and text a lot larger.

With the attribute table in arcgis, as with any program today, I could zoom in by holding the control key and turn the wheel on my mouse. I need to do that in qgis, zoom into the attribute table itself. This is seriously limiting me in using the program, kind of driving me crazy.

Everything I'm looking for talks about how to zoom to a selected feature, I can't seem to word the question right when searching.


texts: Python Classification of multiple NDVI images for difference

texts:
I have several months worth of NDVI images for an area. When I check the values, there seems to be some deviation. I am using python GDAL to attempt to calculate possible areas of concern for inspection.

For Example:

Month     X       Y         V
Jan      55.45  117.25   0.058742
Feb      55.45  117.25   0.0467307
Mar      55.45  117.25   0.153668

And here would be the area

When I compare the NDVI images, some of them seem off. So somehow I need to normalize the datas indexes to find the changes.

I'm kind of stuck and looking for some advice. Perhaps generate VCI's for each image and compare the percentage of change over the area? Each images indexes are correct with itself, but not always relative to the others covering the same area, but I do believe this will work because the min and max values are relative to that image.

Edit:

But the indices are relative only to that image, or so it seems. So the variations within that image are relative to that image only. When I start subtracting to find the changes I end up with some areas that when inspected have not seemed to change.

My objective is to find changes with significant vegetation death and go inspect the cause.


texts: Programmatic Configuration of ArcGIS Server data and layers?

texts: I'm trying to figure out how people in the ArcGIS Server world script configurations.  In the GeoServer world we have the REST Config API, which lets you configure connections to a database or files, publish new layers, update their service metadata, cache the layer, change available and default styles, default projections, etc.  What's the equivalent for ArcGIS Server?  As far as I could tell the REST API there just allows access and manipulation of features and images, but doesn't let you actually publish new data.  Or you can use the admin GUI, but I couldn't find a way to script the types of actions an admin would take.  I've been scouring the docs, so if someone can point me to a reference that would be great.


texts:
Whenever I launch QGIS 2.18.23 all my template print composers launch in the background without any prompting and I'm unable to close them

How can I switch them off or close them?


texts: Batch geocode addresses using a REST API (ArcPy, CSV, Server)

texts:
My stand-alone Python script tries to geocode a CSV file of addresses.
The script successfully geocodes this CSV file when using the *.loc file for this locator but not the published REST API.
ArcGIS Server is federated enterprise server (not portal) version 10.8.1, Python IDLE 3.7.9
Is there something wrong with this call to use the REST API?
arcpy.geocoding.GeocodeAddresses(r"M:\sample.csv",

                             r"https://mycity2.houstontx.gov/pubgis03/rest/services/LOCATORS/COH_ADDRESS_ROAD_LOCATOR/GeocodeServer",
                             "'Address or Place' Address VISIBLE NONE;Address2 <None> VISIBLE NONE;Address3 <None> VISIBLE NONE;Neighborhood <None> VISIBLE NONE;City <None> VISIBLE NONE;County <None> VISIBLE NONE;State <None> VISIBLE NONE;ZIP <None> VISIBLE NONE;ZIP4 <None> VISIBLE NONE;Country <None> VISIBLE NONE",
                             r"M:\BARCtest.gdb\sample5_GeocodeAddresses")

From a broswer, this REST API call works but it uses findAddressCandidates() to get results for just one address.
https://mycity2.houstontx.gov/pubgis03/rest/services/LOCATORS/COH_ADDRESS_ROAD_LOCATOR/GeocodeServer/findAddressCandidates?outFields=*&outSR=4326&f=pjson&SingleLine=1506%20Johnson%20ST
Error tracebacek from the arcpy.geocoding.GeocodeAddresses() call above.
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000005: Could not open the address locator.
Failed to execute (GeocodeAddresses).
Another approach I tried (unsuccessfully) was to use an ArcGIS Server connection file to access the locator.
arcpy.geocoding.GeocodeAddresses(r"M:\sample.csv",

                             r"C:\Users\arcgismanager.HOUTX\AppData\Roaming\Esri\ArcGISPro\Favorites\REST_services_on gisproarc04s_6443.ags/LOCATORS/COH_ADDRESS_ROAD_LOCATOR/GeocodeServer",
                             "'Address or Place' Address VISIBLE NONE;Address2 <None> VISIBLE NONE;Address3 <None> VISIBLE NONE;Neighborhood <None> VISIBLE NONE;City <None> VISIBLE NONE;County <None> VISIBLE NONE;State <None> VISIBLE NONE;ZIP <None> VISIBLE NONE;ZIP4 <None> VISIBLE NONE;Country <None> VISIBLE NONE",
                             r"M:\BARCtest.gdb\sample5_GeocodeAddresses")

This gives the following error.
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Address Locator: Dataset C:\Users\arcgismanager.HOUTX\AppData\Roaming\Esri\ArcGISPro\Favorites\REST_services_on gisproarc04s_6443.ags/LOCATORS/COH_ADDRESS_ROAD_LOCATOR/GeocodeServer does not exist or is not supported
ERROR 000005: Could not open the address locator.
Failed to execute (GeocodeAddresses).
ESRI Technical Support
This ESRI document talks about an ArcGIS Online locator which I'm not using. It also mentions a composite locator. I was using a composite locator but I switched to a non-composite and still got the same error message.


texts:
QGIS 3 is out but doesn't have the plugins I need yet so I'll use QGIS 2.18 a bit longer. I have QGIS 3.2 installed alongside to test the new features.

However, there is now a footer in QGIS 2.18 telling me there is a new version and I should upgrade.

Is there a way to get rid of the yellow footer?


texts: How to subset data for whole county to just portion in particular city using ArcGIS for Desktop?

texts:
I have a shapefile of all of Riverside County. It includes parcels with and without site addresses. I do not have a city limits shapefile. How can I locate the boundary of the city that I need to pull, so I can only grab the attribute data and spatial data for that city and not the WHOLE county?

I tried to do a selection on the site address and used city as a locator, but not all parcels have a site address, which leaves out lots of parcels in that selection.

Can someone provide me with info on how to make a selection on city which will include all parcels, even ones with out a site address? I am using ArcGIS 10.3.


texts: Maintaining spatial integrity of Matlab matrix with missing data

texts:
I am working with MODIS Terra imagery (MOD13Q1).  I am currently trying to compute the autocorrelation structure of the pixels in an MODIS NDVI raster.  The problem is that I have clipped the original MODIS raster to focus on certain regions of interest (ROI).  This means the MODIS raster has lots of missing data pixels in between my ROIs.  I imported the MODIS raster as a netCDF file, then transformed it to a double matrix:

{clear;

%ndvi data
ndvi = double(ncread('j193dzpad.nc', 'NDVI'));
NDVI = ndvi/10000;

%apply no data values
nodata = min(min(NDVI(:,:,1)));
NDVI(NDVI==nodata) = NaN;}

In my case, the minimum pixel value in the MODIS raster is the no data value.

The problem is that when I try to loop through the rows of the MODIS matrix (say to compute my directional autocorrelation values), MATLAB returns NaN.

It's easy to remove the NaN values from the matrix, but doesn't this alter the shape of the matrix?

In my case, I want to maintain the spatial integrity of the data, so I cannot remove the NaN pixels.

Any thoughts on how to work with NaN matrices where one does not want to remove NaN values?

Should I convert my original GeoTiff to another format?


texts: Creating a openstreetmap renderer in android

texts:
I am trying to learn how the openstreetmap is rendered in android. There are already some tools out there that are rendering the maps. Specially I am focusing on offline maps renderer. There is one called osmand.

So, I wanted some suggestions should I go and study how osmand renders the map or do from scratch. I know osmand has done lot more than just rendering. It compresses the maps in obf format. So, I need to study its osmand binary format as well.

So where should I start? I mean should I go for osmand or something else. Osmand have specific implementation.


texts: Numbering adjacent polygons in sequential order

texts: In QGIS is there a way to assign ascending numbers (attribute) to polygons that are located in a row.  The polygons were not created in sequential order.


texts: LEAFLET - geojson properties in bounding box

texts:
Is there a way to show geojson properties which are located in current bounding box?

something like this http://turbo87.github.io/leaflet-sidebar/examples/listing-markers.html
with geojson and its props


texts: Transformation issue with SVG to GeoJSON conversion

texts:
I have a vector image in SVG format. (without georeferences)

I am trying to display it using OpenLayers to display properties of different polygons in it.

I tried converting it to GeoJSON but the resulting render appears to be stretched (ratios are not preserved).

Is there a transformation or a map projection that I am missing or recreating a GeoJSON from scratch the only relevant choice?


texts: Tasseled Cap Transformation coefficient and Bias Value

texts:
I want to perform Tasseled Cap transformation of Landsat 5,7 and Landsat 8 data manually in ArcGIS. Please let me know where from I can get the coefficient values and Bias values for the following formula
TC = Coefficient *DN+Bias
Please share the source too.


texts:
We're migrating from Mapinfo to GIS + PostGIS. I have a file questions about file conversion:

Is there any difference between importing directly from tab to PostGIS, compared with converting tab-shapefile-PostGIS? Do either approaches have accuracy issues? We have hundreds of tab files and there is a mix of lat/long and UTM projections.
What's the best way of converting multiple tab files to PostGIS at once? This is my (poor) attempt at using ogr2ogr:

for %f in (*.tab) ogr2ogr -f "PostgreSQL" PG:"host=myhost user=myloginname dbname=mydbname password=mypassword".

I do have the PostGIS shapefile and DBF loader installed but it only works for shapefiles.


texts:
I am trying to symbolize some polygons with charts, 3 columns each, from the attribute table that has all the fields and rows populated. I have checked for exclusions and queries and there are none. I am only using some simple vbscript for labels but that should not have any effect as it had not in previous maps that I created using the same shapefile only different fields. According to this the problem should be caused by the fields but on the other hand that should also prevent other charts from drawing...

Any ideas?

I am using ArcMap 10.2.2 for desktop.


texts: Removing extra dots when loading Python script

texts:
This is probably a silly question, but I'm having a problem whenever I load a Python script in ArcMap where extra dots are loaded with the script (see below) that I must manually delete prior to running the script. This gets quite annoying, especially when having to run long scripts multiple times. Whenever I write a new script, I just click 'Save As' and specify a Python file, and later load the Python file through the Python window. Is there any way to format the Python files/properties to avoid this from occurring?

... ... # Import system modules
... ... import arcpy
... ... from arcpy import env
... ... from arcpy.sa import *
... ... import math


texts: Convert point to grid: Numbers of cells differ

texts:
In R I am trying to aggregate points data to a grid.
The data for the grid (G-Econ) comes as a csv-file.
In this file each row represent a grid-cell at 1x1 degree resolution and the given coordinates, in longitude and latitude, label the south-west corner of the grid-cell.

In trying to convert these points to a grid I run into the issue that the number of grid-cells in the output (so the eventual grid-raster) does not correspond with the number of cells in the input.

Basically the short version of my question is: How do I solve this?

The simplified G-Econ csv-file can be found here and below I describe step-by-step what I've been doing.

First I load the csv-file, limit the number of observations to only include cells in African countries, and create a spatial object out of the data frame.

# Load data
gecon<-read.csv("gecon.csv",header=TRUE,sep=",",row.names=NULL)
gecon$Cell<-1:nrow(gecon) # Create a grid cell identifier

# Limit data to African countries
require(countrycode)
gecon$cont<-countrycode(gecon$NEWCOUNTRYID,"iso3n","continent",warn=TRUE)
g.africa<-na.omit(gecon[gecon$cont=="Africa",])

# Create spatial object
require(raster)
coordinates(g.africa)<-~LONGITUDE+LAT # Creates spatial object
proj4string(g.africa)=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

summary(g.africa)

From this we see that Africa is covered by 3550 grid cells in the dataset
Next step it to transform the points into a grid:

# Transform points to grid
extent<-extent(g.africa)
extent # Determine number of cols. and rows for grid
length(-26:63) # Number of cols: 90
length(-38:37) # Number of rows: 76
rast<-raster(extent,ncol=90,nrow=76,
         crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
grid<-rasterize(g.africa,rast,g.africa$Cell,fun="last") # Create grid
grid

During this step the number of cells has increased -to 6840- this is due to the inclusion of cells that cover the oceans. I included a cell identifier so I can check how many of the new cells do not belong to a country.

summary(grid)

There are 4053 NAs. Subtract these from the total number of grid-cells and the result is 2787 grid-cells, which is way less than the original number of cells for Africa.
I transform the grid into polygons and this number is confirmed.

africa<-rasterToPolygons(grid) # Create spatial polygons
africa # See nfeatures

So I guess something goes wrong when I rasterize the data, that somehow the interpolation of points to grid-cells results in this reduction.
Due to my inexperience with GIS in R -and in general- I have been unable to solve this issue.

Does anyone have any suggestions on how I could solve this? Or maybe recommendations for a better method?


texts:
I am trying to change geographic reference system images from landsat 8. I created one method to do it in my class. I used the temporary directory as strategy to convert files and clean all row rasters the routine has converted. However I had this error in my script (sh: 1: cannot open TemporaryDirectory: No such file). Does gdal read temporary files? Why does my code show up this error?

def change_src_from_lc8_image(self):

    tmp = tempfile.TemporaryDirectory()
    output_folder_reprojetcted = '{}/{}'.format(self.tmp,self.file_name_targz)

    if not os.path.exists(output_folder_reprojetcted):
        os.mkdir('{}/{}'.format(self.tmp, self.file_name_targz))

    with tarfile.open(self.raster_file_path_targz, "r") as tar:
        tar.extractall('{}/'.format(tmp))

    list_raster_folder = glob('{}{}'.format(tmp, '/*/*TIF'))

    for tif in list_raster_folder:
        command = "gdalwarp {img_src} {img_output} -s_srs EPSG:32624 t_srs EPSG:4674".format(img_src=tif, img_output=output_folder_reprojetcted)
        os.system(command)

    tmp.cleanup()


texts:
I created a project in QGIS 2 and put in an Atlas. I then realised that QGIS 2 only allows 500 objects in the Atlas preview, which is no good to me (I've got about 2000+) so copied the project, opened it in QGIS 3 and got it working.

I was having no problems with QGIS 3 until the last step, which meant I had to create 9 extra print composers (or "layouts" in QGIS 3). These are essentially copies of the first one but with minor changes (logos, wording etc), so 10 in total.

Now, whenever I try to open the project in EITHER version of QGIS (since I also went back and created them in the QGIS 2 version in case anybody needs that), it takes forever. In QGIS 2 you can see it going through the print composers and presumably loading them into memory, and I assume the same is happening with QGIS 3 - only I thought a feature of QGIS 3 (I'm on 3.4.5) was that print composers / layouts weren't meant to be loaded on startup.

Does anybody know if there's a way to stop them loading on startup? I still need them so deleting isn't an option.


texts: How to retrieve metadata from a feature layer in ArcGIS Server 10 using Flex

texts: I need to retrieve the metadata information from a feature layer. I have the web services available and I can run queries to retrieve the feature attributes, but I don't know any method to retrieve the metadata. Is there any sample or code snippet that I can use to understand the method?


texts: Understanding difference between Amazon Web Services and GeoServer?

texts:
I'm learning how to develop a web mapping application for multiple GIS layers (vector and raster). As I'm learning about each component when configuring every element works (GIS desktop, RDBMS, web server, geospatial server, tiling, APIs), I'm hitting a roadblock in conceptually understanding how to choose and implement a server system.

AWS comes highly recommended from colleagues, but I'm getting confused how it works. I know it's a cloud server, so does this mean I can store/host my GIS data (tiles, WMS) on there? How is it different from GeoServer?

I can't get past understanding the functionality of those two systems.


texts: Difference in raster size using raster2pgsql

texts:
I use the PostGIS-build-in raster2pgsql to load SAGA-raster into the db. I tried to grid them into 100x100-tiles, just as seen in the code:

  /usr/lib/postgresql/9.2/bin/raster2pgsql -s 32632 -a -F -Y ~/Desktop/examplerast.sdat -t 100x100  import.exampletable | psql testrun -h localhost -p 5432 -U dbmaster

The strange thing is, that tiles in PostgreSQL/PostGIS have doubled dimensions: 200x200; as can be seen here:

  select st_asewkt(st_envelope(rast)) from grids.exampletable limit 1;

  SRID=32632;POLYGON((589203.15 5954785.35,589203.15 5954787.35,589205.15 5954787.35,589205.15 5954785.35,589203.15 5954785.35))
  (1 row)

In order to achieve the desired 100x100-tiles I have to load them with a dimension -t 50x50

Why this doubling of dimensions? How to get rid of it? Thanks in advance!


texts: Add Field doesn't work on ArcPy (ArcGIS) in a loop

texts:
I am trying to add a field to each shapefile inside several folders in ArcGIS with Arcpy. I am executing the following code.
import arcpy, os
# Overwrite pre-existing files
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
regions = ["China","India","Euro","ArabPersia"]
years = ["622","800","900","1000","1100","1200","1300","1400","1500","1600"]
for region in regions:
    for year in years:
        folder = r'C:/Users/Dropbox/Trade_Diaspora/LCP/LCP_gvtyCntrs'+"//"+year+"//"+region+"//"+'polylines'
        for root, dirs, files in os.walk(folder):
            for name in files:
                if name.endswith('.shp'):
                    shpName = os.path.splitext(name)[0]
                    absFile = os.path.abspath(os.path.join(root,name))
                    arcpy.DeleteField_management(absFile,"cityid")
                    arcpy.AddField_management(absFile,"cityid","SHORT")
                    arcpy.CalculateField_management(absFile, "cityid",shpName,"VB")
                    print("complete " + region + " " + year + " "+ shpName + " field addition")
        print("Completed All  "+region+ " " +year)

The problem is that it works sometimes for some random times. In some cases, it loops through like 50-60 files before throwing an error. Sometimes it throws an error after 2 loops.
What I have tried:

Tried to run for a single folder (No first two loops)
Tried with adding empty values to options in the AddField command.
Tried to run both on ArcMap's python as well as a stand-alone script.
Started without DeleteField command for files that do not have the field for sure.
Tried with just AddField not the CalculateField command.
Tried with adding another field (also changed names)
Switched my internet connection to see if there is any issue (ArcGIS needs to run on VPN)

None of these worked in my case. I am using ArcGIS 10.6, Python 2.7.14, 32 bit (can't use 64 bit because I need to convert these into Excels later), Windows 10, with 32GB RAM.
Every time, I get the same error in the AddField command:
Traceback (most recent call last):
  File "C:\Users\Dropbox\Trade_Diaspora\LCP\LCP_gvtyCntrs\addfieldonly.py", line 16, in <module>
    arcpy.AddField_management(absFile,"cityid","SHORT")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcPy\arcpy\management.py", line 3424, in AddField
    raise e
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (AddField).


texts: How to add a polygon within another polygon as opposed to on top of it?

texts:
Sorry for the confusing question. I have a shapefile and I need to draw a polygon in the middle of it. I know I can use the "create features" function, but that creates a polygon on top of the first one, whereas I need it to cut into the main polygon and reduce the area of it, rather than being an independent feature sharing the same space. The same way that the "cut polygons" function divides up the original polygon, but I can't always start from an edge.

Does that make sense? It seems like such a simple procedure but I can't find any way to do it.

I've attached an image to help explain. When I select the main polygon, it doesn't highlight the edges of the smaller polygons within it, meaning they're not altering the shape/area of it.


texts: In ArcMap 10.6 how can I turn off the 'snapping' associated with drawing rather than features?

texts:
I have an Arcmap MXD where I have created annotation and am adjusting its placement, but it seems to keep snapping to an unseen grid and I can't place things exactly where I want them.
I am looking for a way to turn off this behaviour so the labels stay where I drop them.
I have done some internet searching but keep coming up with answers related to feature snapping to other features, or to the guides at the edge of the neatline, whereas I do not have those guides in my MXD and am not snapping to another feature, but rather whatever invisible grid is bring applied to the MXD.  Any help?


texts:
I have installed GeoServer 2.16 on a brand new Windows 10 laptop at C:\Program Files\geoserver-2.16.0.  I can open the interface on http://localhost:8080/geoserver and am logged in as admin.  I am working through the Getting Started section of the user manual, and have reached 3.2 Publishing a shapefile.

I have unzipped a folder of shapefiles into the data directory at C:\Program Files\geoserver-2.16.0\data_dir\data\MyShapefiles.  I have created a workspace.  When I try to create a Store, I choose Shapefile, get the New Vector Data Source page, select the workspace I created, and enter a name.

When I try to choose the Shapefile location, I click Browse but I don't see my directory of shapefiles.  I am in the Data directory, and it displays gwc/, logs/, security/, styles/, and workspaces/.  I can't find my directory nor any of the sample data that comes with GeoServer installation.


texts:
I've encountered some performance issues when developing in ArcPad 10. The first is the loading of the map file. It takes 1 minute with sample data, and once opened each form takes 1 to 2 sec to be displayed, especially when the form contains comboboxes filled from other tables. But it takes 5 sec for some complex geometric calculations including search in all the features. (I think that this delay could increase with data size.)

I think that this performance is due to a large geodatabase including 20 layers with just 6 visible layers. The invisible layers are more like related tables. I used related tables in this way to avoid some design limitations in exporting related tables to ArcPad.

The Question:
Can a large number of invisible layers reduce ArcPad performance?
What could I do to improve the performance of this ArcPad Application?


texts: Testing if Geometry is convex using JTS?

texts:
I'm using the Java Topology Suite and I'd like to efficiently check if a Geometry (a Polygon in fact) is convex or not.

I currently do if (geometry.convexHull().equals(geometry)) but obivously that is ... suboptimal :)

I could implement it myself with the early-interruption condition on the gift-wrapping algorithm, but I can't imagine that this isn't implemented in JTS already.


texts: Settings for GeoServer: large WMS png request failing

texts:
We recently upgraded from GeoServer 2.1 to the 2.5-SNAPSHOT, both from official OsGeo Centos packages, git revision, c7886eb21885c88d40029432c3178f88931ed820

We are getting the dreaded:

  Rendering process failed
  java.lang.OutOfMemoryError: Java heap space
  Java heap space

errors on processing. The WMS png image we are requesting is 2600 x 2600 pixels and worked on the older version of GeoServer. If we reduce the width/height below 2000 x 2000 pixels, the image is returned. This would seem to imply some cutoff at  around 15Mb, assuming 3/4 bytes per pixel, and 4 million pixels at 2000 x 2000.

We have tried setting the config variable <maxRequestMemory>50000</maxRequestMemory> in web.xml in the Geoserver config directory as suggested in the docs.

We have also tried setting Xmx via JAVA_OPTS="${JAVA_OPTS} -XX:MaxPermSize=128m -Xmx768m" in /etc/tomcat6/tomcat6.conf as suggested in the docs.

Running,

ps ax|grep tomcat

yields

  /usr/lib/jvm/jre/bin/java -XX:MaxPermSize=128m -Xmx768m -Djavax.sql.DataSource.Factory=org.apache.commons.dbcp.BasicDataSourceFactory -Djava.awt.headless=true -Xms256m -Xmx768m -Xrs -XX:PerfDataSamplingInterval=500 -XX:MaxPermSize=256m -Dorg.geotools.referencing.forceXY=true -DGEOEXPLORER_DATA=/var/lib/opengeo/geoexplorer -XX:MaxPermSize=128m -Xmx768m -Djavax.sql.DataSource.Factory=org.apache.commons.dbcp.BasicDataSourceFactory -Djava.awt.headless=true -Xms256m -Xmx768m -Xrs -XX:PerfDataSamplingInterval=500 -XX:MaxPermSize=256m -Dorg.geotools.referencing.forceXY=true -DGEOEXPLORER_DATA=/var/lib/opengeo/geoexplorer -classpath :/usr/share/tomcat6/bin/bootstrap.jar:/usr/share/tomcat6/bin/tomcat-juli.jar:/usr/share/java/commons-daemon.jar -Dcatalina.base=/usr/share/tomcat6 -Dcatalina.home=/usr/share/tomcat6 -Djava.endorsed.dirs= -Djava.io.tmpdir=/var/cache/tomcat6/temp -Djava.util.logging.config.file=/usr/share/tomcat6/conf/logging.properties -Djava.util.logging.manager=org.apache.juli.ClassLoaderLogManager org.apache.catalina.startup.Bootstrap start

which shows 768m for Xmx, which should be vastly in excess of the requirements for a 2600x2600 pixel png. This is against a dev server with no other users or requests against it. Clearly we are missing something.


texts: smallest number of cells classified in a raster

texts:
I am using ArcGIS 10.3 and I ran con statement in the ModelBuilder within ArcGIS. I have an output raster (2m*2m) classified into four different classes. The raster has small pockets of 4 classes. I want to know what is the smallest number of cells that ModelBuilder has been able to categorize into a particular class.

Is there any specific command to figure this out or randomly we can say it can pick even a single cell that falls into a separate class?


texts: OpenLayers custom tileGrid for scanned map image - How to calculate resolutions?

texts:
I'm using OpenLayers 3 and attempting to display a layer based on a custom tile grid in a custom projection (EPSG:4202) on a standard projection map (EPSG:3785).  The layer is essentially a photo of a map that I took, and then split it into a tile grid of images.  I have successfully created the directory structure and image files, and OpenLayers can see them and display them, putting all the tiles together in the right sequence relative to each other.
However, the layer is displayed much, much smaller than it should be.  The origin is in the right place, but the image is tiny... all shrunk up towards the origin.
I guess I'm calculating the resolutions incorrectly.  This also implied by the insane number of tiles that OpenLayers wants to fetch (and when zoomed out, the OpenLayers error that appears in the JavaScript console when zoomed out that says:  Assertion Failed: reasonable number of tiles is required).
Eventually it displays a small rectangular image  where the origin (top-left) of this image is in the correct location.  But the image should be much larger and should extent much further down and to the right.
Complete code is below.
<!doctype html>

    <link rel="stylesheet" href="ol/ol.css" type="text/css" />
    <link rel="stylesheet" href="map.css" type="text/css" />

Here's how I calculated the resolutions (this is pseudo-code, not JavaScript):
resolution = ( imageExtentInLayerProjection[2] - imageExtentInLayerProjection[0] ) / imageWidth
Each subsequent resolution is:
resolution /= 2
(imageExtentInLayerProjection is the extent of the intended layer's spatial reference system that the original map image should cover.  It is an array of the form:  [minX, minY, maxX, maxY].  The values were taken from the map graticule labels.  The actual values are:  [146.5, -41.5, 147, -41]).
What am I doing wrong here?


texts: Individual tree segmentation using R-package lidR

texts:
I have some questions about individual tree segmentation in lidR. I'm looking at using the function lastrees with algorithm watershed.

My first question is how to apply it for large areas. How big areas or how big datasets (number of points) can be handled in one call to lastrees?

My second question is, if there is an upper limit on the data size, is there then a way to handle bigger areas by combining the resulting segments? Is that included in lidR or should the user implement this?

Finally, is there any publication where this has been used and validated for segmentation and forest estimation and does it include suitable parameter settings and ways to define a chm for the watershed?


texts:
I noticed that PyQGis has QgsMapToolPan() and QgsMapToolZoom() functions, and the cursor shape on the canvas will change to a hand (for pan), magnifying-glass(+/-) (for zoomIn / zoomOut).

I have a toolbar for Select (arrow icon), Marquee (bounding box icon), besides the Pan, ZoomIn and ZoomOut tool buttons.

Currently I created the Select & Marquee tools as QgsMapTool(canvas) and setAction with the appropriately ui widget buttons.

How can I change the cursor shape to arrow when the Select tool is activated, and cross when the Marquee tool is activated?

(I noted the default cursor shape for QgsMapTool() is Cross icon which I could use for the Marquee tool.)


texts: Symbology for Open Street Map [OSM] Data in MapInfo [tab] format

texts:
I have shapefiles for the different OSM layers which I have translated into TABs but I am now trying to symbolise them in a similar style to OSM online.

Is there a way to easily apply the OSM symbology to MapInfo TAB files?


texts: Recalculate XY Extent Option Not Available

texts: I am unable to recalculate the XY extent of a feature class. I am attempting to do so within the Feature Extent tab within the Properties window accessed from ArcCatalog. The Recalculate button, however, is greyed out. The feature class is within an ArcSDE geodatabase and I am not the DBO. Any ideas on why this option is not available?


texts: Splitting text attribute using conditional function QGIS 2.18.13

texts:
i want to split a text column with city names by selecting only the names like 'world (planet)' and not like 'world planet'. Tried before with left( "Name", strpos( "Name" ,' ')) update existing field and results are fine only for the names 'world (planet)' but the other names like 'world planet' disappeared from column, So tried with conditional function case when 'names'like '%(%' then left( "Name", strpos( "Name" ,'')) end (if this its a correct form of expression).Is there a way to split the names which is included the '()'?
(im new user in QGIS.)


texts: ST_Intersects incorrectly returning false in PostGIS

texts:
I have a shapefile created from ArcMap that contains about 400 geometry's. I've added this shapefile into PostGIS using the shp2pgsql tool. I am trying to use the ST_Intersects function to determine if a given point, using Lat and Long, is within one of the 400 geometry's.
I have a point that I know to be within one of the shapes, but ST_Intersects is returning false for this. The record where cocnum equals "MA-500" contains the point Lat: 42.3581, Long: -71.0636.

                    SELECT
                        cocname,
                        cocnum,
                        geom,
                        st_intersects(geom, point(-71.0636,42.3581)::geometry)
                    FROM
                        public.coc
                    WHERE cocnum like 'MA-500'

However, I get the following result:

When I create a KML of the geometry using st_asKml(geom) and then insert that into google maps, along with the point google.maps.LatLng(42.3581, -71.0636) I get the following result:

Any ideas on why ST_Intersects is returning false for this?


texts: Viewing the data without any background imagery in Cesium

texts: Is there a way to completely disable the background imagery in cesium? I would like to see only my data without any background imagery. Is this possible?


texts: Getting TopologyException: Input geom 1 is invalid which is due to self-intersection in R?

texts:
The 'TopologyException: Input geom 1 is invalid' self-intersection error which arises from invalid polygon geometries has been widely discussed. However, I haven't found a convenient solution on the web that solely relies on R functionality.

For instance, I have managed to create a 'SpatialPolygons' object from the output of map("state", ...) following Josh O'Brien's nice answer here.

library(maps)
library(maptools)

map_states = map("state", fill = TRUE, plot = FALSE)

IDs = sapply(strsplit(map_states$names, ":"), "[[", 1)
spydf_states = map2SpatialPolygons(map_states, IDs = IDs, proj4string = CRS("+init=epsg:4326"))

plot(spydf_states)

The problem with this widely applied dataset is now that self-intersection occurs at the point given below.

rgeos::gIsValid(spydf_states)
[1] FALSE
Warning message:
In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point -122.22023214285259 38.060546477866055

Unfortunately, this problem prevents any further use of 'spydf_states', e.g. when calling rgeos::gIntersection. How can I solve this issue from within R?


texts:
I'm using QGIS 2.18.9.
Using the Handle Bad Layers Window, I find that I can't change the window size or the size of the individual columns.  Sometimes I have a long path and/or file name and see only the first few characters in the relevant column.  I might have several layer files with the same first few characters and hence don't know what file to browse for.
Is there some way of enlarging the window and columns?
I have just discovered that if I double-click on the visible text in the Datasource column, the complete filename is displayed.


texts: Unable to update features in tables without FID: How to use gid field as the FID column with ogr/gdal python?

texts:
I'm trying to update a PostGIS table using the SetField() method in python. The table is imported from a shapefile  with pgsql2shp so it has a gid field, which can be used as FID or primary key.

However, when executing SetField(), I get an error saying:

Error 1: Unable to update features in tables without a recognised FID column.

I tried different things. The only thing that allows updating is to externally add and populate a new column called ogc_fid. Even adding a 'fid' column won't do.

Is there a way to tell ogr/gdal (as of gdal2.1.3) in the python code itself that an existing field such as gid is the FID column?


texts: Identify name of archive class using arcpy

texts:
I am trying to use arcpy to identify the name of the archive class for a number of feature classes, but cannot find a function that will output the name of the archive class, only whether the feature class has archiving enabled.

import arcpy

gdbConn = r"Database Connections\ScriptTest.sde"
arcpy.env.workspace = gdbConn

queryResult = ["List of all my feature classes"]

for row in queryResult:
    gisLayer = row
    desc = arcpy.Describe(gisLayer)
    if desc.isArchived:
        print "{} has Archiving Enabled".format(gisLayer)

        gisLayer_Archive = "" # Identify Archive Class here

        arcpy.FeatureClassToFeatureClass_conversion(gisLayer_Archive, r'D:\temp\TempOutput.gdb', gisLayer_Archive)

My intention is to take a copy of each archive class and save them in a file geodatabase.

Is there a way to find the name of the archive class for archived feature classes in a geodatabase so that I can copy those archives out?

I could possibly assume that all archive classes have a _H suffix, however it is possible that some of these have changed over time.

ArcGIS for Desktop/Server 10.3.1
MS SQL Server 2012


texts: How to manually break statement execution in Python window of ArcMap?

texts:
How does one break the execution of a long-running process in the python console of ArcMap?

For example, how do you break a long iteration like the one below once it is started?

for i in range(1, 50):
  # some long process which takes over a minute to complete

I want to stop the process at any time I want by pressing a keyboard combination, clicking a button or something similar while the code I have input in the console window is running.


texts: Merge files starting with same characters using ArcPy

texts:
I am trying to merge files starting with same 2 characters as below;

11111  - merge (output name : 11.shp)
11112  - merge (output name : 11.shp)
11113  - merge (output name : 11.shp)
22112  - merge (output name : 22.shp)
22153  - merge (output name : 22.shp)

But my python script' result merge together if "11" is including in the filenames(like 4)).
I only want 11***.shp(1)~3))   but A Script merges 11.shp (merge 1~4 as 11.shp)
========================================================================
import sys, os, arcpy, fnmatch
reload(sys)

arcpy.env.workspace  = r"F:\a"

workspace = arcpy.env.workspace
outdir = r"F:\b"

List = []

for dirpath, dirnames, filenames in os.walk(workspace, topdown=True):
    for filename in filenames:
        List.append(os.path.join(dirpath, filename))
    break

patterns = set([os.path.basename(fl)[:2] for fl in List])

for pattern in patterns:
    group = [fl for fl in List if fnmatch.fnmatchcase(fl, '*' + pattern + '*.shp')]
    output = os.path.join(outdir, pattern)
    arcpy.Merge_management(group, output)

========================================================================
How can I fix it?


texts: Calculate areas in ha (hectares) with gdal?

texts:
I would like to calculate using gdal certain areas in ha (hectares) for a reclassified raster. Can somebody help me with this?
Here is the reclassification:

#RBR reclassification to USGS burn severity values

arr = np.array(RBR)

level1 = (arr >= -0.500) & (arr <= -0.251) #Enhanced Regrowth, high (post fire)
level2 = (arr >= -0.250) & (arr <= -0.101) #Enhanced Regrowth, low (post-fire)
level3 = (arr >= -0.100) & (arr <= 0.99) #Unburned
level4 = (arr >= 0.100) & (arr <= 0.269) #Low Severity
level5 = (arr >= 0.270) & (arr <= 0.439) #Moderate-low Severity
level6 = (arr >= 0.440) & (arr <= 0.659) #Moderate-high Severity
level7 = (arr >= 0.660) & (arr <= 1.300) #High Severity

levels = [level2, level3, level4, level5, level6, level7]
reclass_values = [-1, 0, 1, 2, 3, 4]

reclassed_arr = np.where(level1,-2, arr)
for level, val in zip(levels, reclass_values):
    reclassed_arr = np.where(level, val, reclassed_arr)

And this is the burn area:

#Calculate burn area

reclass = reclassed_arr
k = ['Unburned hectares', 'Low severity hectares', 'Moderate-low severity hectares', 'Moderate-high severity hectares', 'High severity']
for level, val in zip(levels, reclassed_values):
    x = reclassed_arr[reclassed_arr == level]
    l = x.size*0.4

I got this error:

----> 8 print("%s: %.2f" % (k[reclassed_arr],1))
TypeError: only integer scalar arrays can be converted to a scalar index

I would like to do something like this code to obtain the size of the different category areas.

reclass = reclassify(clipped_ds_rasterband.ReadAsArray())
k = ['Unburned hectares', 'Low severity hectares', 'Moderate-low severity hectares', 'Moderate-high severity hectares', 'High severity']
for i in range(1,6):
  x = reclass[reclass == i]
  l= x.size*0.04
  print("%s: %.2f" % (k[i-1], l))

At this point I get the number of pixels for each category.


texts: Finding upstream and downstream distance between points along river network polyline

texts:
I have a network of polylines making up a stream network (from the National Hydrology Dataset). I also have 2 point files: one of areas of interest and another of locations of dams within the network. I'd like to find a way to calculate 1) the distance to the closest upstream dam for each of the areas of interest and 2) the distance to the closest downstream dam.

I would prefer to find a way to do this in ArcGIS but have also been looking into using R or Python script.


texts:
I have noticed with two mobile devices (Samsung S2+ and OnePlus 5) that the GNSS sensor picks up QZSS satellites (PRN 193 and 194 [only for OnePlus 5]), even while using these devices over Northern parts of India. Strangely, I have not picked up these satellites with some dedicated GPS receivers (Garmin GPS12, Garmin GPSMap 78s and Magellan GPS315) and even some mobile devices (Motorola Defy+, Motorola Moto G3, Nexus 7C, Nexus 7 2013), where as the QZSS signal is supposed to be compatible with traditional GPS.

This webpage (official, I think) shows a graphical depiction of Main area where QZSS is avalable as:

So, is/are there any utility of my device picking up QZSS satellites even while being in Northern parts of India?


texts:
I am trying to run saga tool however I always received errors.

I am using Ubuntu Xenial Xerus (16.04), Qgis 2.14.6 and Saga 2.2.7.

Log (OGR) error:
2016-09-07T15:22:28 1   Data source is invalid ()

See images of errors:

Log (Processing) Error:
Error: input file [/tmp/processing25e6c401b76f41efb6bb44365b964fdf/4ee2c30104ca4cf39033d19cdb51167e/srtmriachosantanaSIRGAS200023S.sgrd]
            Error: DEM
            ____________________________

            #####   ##   #####    ##
            ###     ###  ##       ###
            ###   # ## ##  #### # ##
            ### ##### ##    # #####
            ##### #   ##  ##### #   ##
            ____________________________

            SAGA Version: 2.2.7

            ____________________________
            library path: /usr/lib/x86_64-linux-gnu/saga/
            library name: libta_preprocessor
            library     : Preprocessing
            tool        : Fill Sinks XXL (Wang & Liu)
            author      : Copyrights (c) 2007 by Volker Wichmann
            processors  : 8 [8]
            ____________________________

            Load grid: /tmp/processing25e6c401b76f41efb6bb44365b964fdf/4ee2c30104ca4cf39033d19cdb51167e/srtmriachosantanaSIRGAS200023S.sgrd...
            failed

            Usage: saga_cmd ta_preprocessor 5 [-ELEV <str>] [-FILLED <str>] [-MINSLOPE <double>]
            -ELEV:<str>         DEM
            Grid (input)
            -FILLED:<str>       Filled DEM
            Grid (output)
            -MINSLOPE:<double>  Minimum Slope [Degree]
            Floating point
            Minimum: 0.000000
            Default: 0.100000
2016-09-07T15:22:28 2   Error loading result layer:
            Traceback (most recent call last):
              File "/usr/share/qgis/python/plugins/processing/gui/Postprocessing.py", line 76, in handleAlgorithmResults
                out.name))
              File "/usr/share/qgis/python/plugins/processing/tools/dataobjects.py", line 199, in load
                + '\nCheck the processing framework log to look for errors')
            RuntimeError: Could not load layer: /tmp/processing25e6c401b76f41efb6bb44365b964fdf/b5e9bd717ed34b15a08ebaf91e8bee8b/FILLED.sdat
            Check the processing framework log to look for errors


texts: Free/Open-Source Archaeological Geophysics Program?

texts: Does anyone know a free/open source software or a GIS plugin for working with archaeological geophysics data, especially for data from a magnetic survey (magnetic gradiometer map)?


texts: How to style markers so they dynamically offset if overlapping at different zoom levels?

texts:
I have created a map using Carto. The map has bespoke images as icon markers. I would like to style the map so that, if the icons overlap,
  they both render in an offset position so that the full icons are
  visible, perhaps with a line to the actual co-ordinates. Ideally a
  function for this exists and can handle any overlap dynamically, if
  not some manual tips on offset styling will help.

  An example of the issue I have is here: https://www.mapbox.com/mapbox-gl-js/example/custom-marker-icons/ If
  you zoom out far enough on the map, the icons will begin to overlap,
  which is the behavior I wish to prevent.

  Clustering is not an option as I have unique images that need to be visible.

  Spidering/Spiderfy (http://jawj.github.io/OverlappingMarkerSpiderfier-Leaflet/demo.html)
  looks close to what I am after, but this only splits the icons on
  click rather than as standard.

Moving on from the original question above, i have managed to create hard-coded offsets for overlapping features with the below sql:

    SELECT cartodb_id, tree, capacity, the_geom_webmercator FROM tablea  where tree NOT IN ('Oak','Willow')
UNION SELECT cartodb_id, tree, capacity, st_translate(the_geom_webmercator,20000,0) as the_geom_webmercator FROM tablea where tree = 'Oak'
UNION SELECT cartodb_id, tree, capacity, st_translate(the_geom_webmercator,10000,10000) as the_geom_webmercator FROM tablea where tree = 'Willow'

This brings in all other elements in tablea and then styles specific elements with an offset on the X,Y bringing them all together as a union select.

Unfortunately this creates a problem. The hard coded offset is fine at certain zoom levels but looks odd in others. Ideally I want a zoom specific sql query. From what I can tell this doesn't exist in Carto, and you can only style the CSS at specific zoom levels?

Mapbox have created a postGIS function that runs the below:

( SELECT * FROM roads
  WHERE Z(!scale_denominator!) >= 12
) AS data

Does anyone know of the Carto equivalent to this? Or perhaps advice on an alternative method?


texts:
I am trying to request features of a wfs layer and then zoom to them. The wfs layer is defined as follows:

wfsLayer = new OpenLayers.Layer.Vector("University", {
      strategies : [new OpenLayers.Strategy.Fixed(), cluster_strategy],
             styleMap : new OpenLayers.StyleMap({
                 "default" : gebaeudeNrStyle,
                 "select" : gebaeudeNrStyle
             }),
      protocol : new OpenLayers.Protocol.WFS({
                 url : "/geoserver/wfs?",
                 featureType : "Gebaeude",
                 featureNS : "namespace",
                 geometryName : "the_geom",
                 srsName : "EPSG:4326",
                 version : "1.1.0"
      })
});

Here the selectFeature function, which requests features of the above wfs layer by an id and zooms to it.

function selectFeature(id) {

    var feature = wfsLayer.getFeaturesByAttribute("Gebaeude_Nr", 03);

    if (feature) {
         map.zoomToExtent(feature.geometry.getBounds());
    }
}

In this example, i'd like to zoom on the third id (Gebaeude_Nr), which inexplicably does not work. getFeaturesByAttribute() delivers an empty array. The feature.length is zero and firebug returns an error "TypeError: feature.geometry is undefined".
However, the underlying database table for the featureType "Gebaeude" contains a column "Gebaeude_Nr" and an entry "03".

How can i search for a WFS feature and zoom to it in OpenLayers?
Thanks in advance

Side note:

for (var i = 0; i < wfsLayer.features.length; i++){
    alert(wfsLayer.features[i].geometry);
}

works, which gives me outputs like POINT(895866.7304909064 6849316.098470407) for every feature. However, i haven't found out yet how to get all column names of the layer in order to compare the results with the database. This would allow to check if the above column "Gebaeude_Nr" exists for getFeaturesByAttribute("Gebaeude_Nr",03);


texts: Map walking distance/time from transit stops

texts:
I'm working in QGIS and would like to show how far from a bus stop a pedestrian could walk along the street network in any direction for a specified set of walking times (5,10,20 mins etc).

I have my street network as a modified OS ITN Road Centre Line layer and bus stop locations as points.

Here's what i'm aiming for: http://www.humantransit.org/2010/05/culdesac-hell-and-the-radius-of-demand.html rather than isochrone/heat maps.

Many thanks


texts: Trying to get major roads using Name of the road

texts:
I am trying to use OverPassApi to get a road Geopoints of Major/Minor roads in USA.

Here is the query

<osm-script output="json">
  <id-query {{nominatimArea:USA}} into="area"/>
  <query type="way">
  <has-kv k="highway" regv="motorway|trunk|primary|motorway_link|trunk_link|primary_link"/>
<has-kv k="name" regv="West Side Highway"/>
      <area-query from="area"/>
  </query>
  <union>
    <item />
      <recurse type="way-node"/>
  </union>
  <print mode="body" order="quadtile"/>
</osm-script>

I am successfully trying it on http://overpass-turbo.eu. With few exceptions, in general this query is ok. If anyone can make it better it would be helpful.

My question is, I want to convert it to OverPass QL Using

http://overpass-api.de/query_form.html

It gives error at

<id-query {{nominatimArea:USA}} into="area"/>

How can i solve this. Secondly, If i remove this line, and try to run the same coverted code QL on query

NSString *urlString = @"http://overpass-api.de/api/interpreter?data=[out:json];way[\"highway\"~\"motorway|trunk|primary|motorway_link|trunk_link|primary_link\"][\"name\"~\"%@\"](area.area);(._;node(w););out body qt;";

%@ = West Side Highway

it is giving me Empty Array.

What is wrong?


texts: Multicriteria Analysis for flood risk analysis

texts:
I have used a multi-criteria analysis which uses factor weighting for flood risk identification in a river basin using ArcGIS 10.1. My factors are rainfall, drainage, density, slope, and land use.

But the result of the analysis shows non-flooded area as flooded and flooded areas as not flooded. (I have a digitized basemap which shows the flooded area of that river.) What should I do to correct it and have results same as the flooded area base maps published by our countries disaster management center?

Please tell me how to have at least 90 percent accuracy of my result when compared with the base maps which show the flooded areas.


texts:
I have a raster whose Layer properties is as follows:

    Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
D:/STUDY/MSC/Thesis/GIS works/temp/LU_New_region.tif
AREA_OR_POINT=Area
Authors=Sophie Bontemps ,Eric Van Bogaert ,Pierre Defourny
Copyright=Copyright Â©UCL Geomatics, BELGIUM 1999-2010
Generate by=gdal_mean
process begin time=2010-12-23T09:49:37
process files=CL5_GLOBCOVER-L5_CLASSIF_2009_V2.3.20101220.tif
process finish time=2010-12-23T09:57:38
Band 1
STATISTICS_MAXIMUM=220
STATISTICS_MEAN=111.55999946952
STATISTICS_MINIMUM=11
STATISTICS_STDDEV=76.782603856836
Dimensions
X: 5589 Y: 2967 Bands: 1
Origin
76.3157,34.5213
Pixel Size
0.00277774,-0.00277774
No Data Value
0
Data Type
Byte - Eight bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=45 +datum=WGS84 +units=m +no_defs
Layer Extent (layer original source projection)
76.3156519612749094,26.2797412081658877 : 91.8404507576846925,34.5213010630500889
Band
Band 1
Band No
1
No Stats
No stats collected yet

I have another raster whose Layer spatial reference system is as follows:

 Driver
GDAL provider
GTiff
GeoTIFF
Dataset Description
D:/STUDY/MSC/Thesis/GIS works/GAV_X_combined_final.tif
AREA_OR_POINT=Area
Band 1
STATISTICS_MAXIMUM=90
STATISTICS_MEAN=37.946288072175
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=29.945215505785
Dimensions
X: 2386 Y: 4079 Bands: 1
Origin
251326,3.20264e+06
Pixel Size
30.3773,-30.3786
No Data Value
-3.40282e+38
Data Type
Float32 - Thirty two bit floating point
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=45 +datum=WGS84 +units=m +no_defs
Layer Extent (layer original source projection)
251326.4686199999996461,3078724.4376099999062717 : 323806.8145299999741837,3202638.6342000002041459
Band
Band 1
Band No
1
No Stats
No stats collected yet

They don't superimpose. On the fly projection is OFF.
I can see the layer extent is not near in both of them. Is it causing the problem ?


texts:
I deleted SiteAddressPoint in the enterprise geodatabase.  Now when I try to replace it I get an error that says it already exists.  I have disconnected all connections and cleared locks. I have gone through all the REST services and stopped them and made sure "Lock schema" was unchecked.  Here is what I have tried and the result:
Creating a new empty feature class.  It will not let me proceed and says the name already exists.
Importing Feature class to feature class Fails and says the name already exists and It cannot find the feature class. Error 000258  Overwrite in geoprocessing is checked.
Python with arcpy.DisconnectUser(sde, 'ALL') and arcpy.ClearWorkspaceCache_management()  I get ERROR 000601: Cannot delete xxx.DBO.SiteAddressPoint.  May be locked by another application. Failed to execute (CopyFeatures).
I have tried in ArcMap 10.5.1, ArcCatalog, ArcGIS Pro 2.7 and the Python script both in ArcMap and in IDLE.
It is ArcGIS Server 10.5.1 and SQL Server.
How can I make a new SiteAddressPoint in the database?


texts: Altitude information for European highways

texts:
I need to calculate the gradient in major European highways. Fot this I need altitude data for each highway, let's say every 100 meters.
Are there such kind of data publicly available? If not, I am also interested in commercial ones.

I have already tried to use DEM GeoTIFF images, but since the DEM depicts the "bare earth", they are not accurate enough for my purpose.


texts: Using Wildcards with ArcMap Data Driven Pages

texts:
Does anyone know if it is possible to achieve the effects of a wildcard match using data driven pages (DDP) in ArcMap?

For one of the tables that I want the DDP index value to match, the target field has a comma-separated list of values. As long as the DDP index value is in that list, I would like the table to filter in/out those rows.

For example, if the DDP index value is 'Y' and the corresponding field value in the other table is 'X,Y,Z', I would like that row to evaluated as matched. I've tried setting the DDP index values to include wildcard characters (*, %) but that has no affect.

The intended purpose of this is for real-time interaction in ArcMap, so I don't think a solution using arcpy would be helpful. However, if there is some other programming option that is relatively simple, I would be certainly be interested in that.


texts: How to pull in and push out text browser?

texts:
I am using Qt designer to create a plugin in Qgis3. I want to add description of a tool beside the tool. For this I have tried using 'Text Browser' by connecting it with button. I am able to hide and display the description(A Text Browser) using following code:

    def hide(self):
    if self.dlg.textBrowser.isHidden()==True:
        self.dlg.textBrowser.setVisible(True)
    else:
        self.dlg.textBrowser.hide()
    #connect the action using Push button
    self.dlg.pushButton_5.clicked.connect(self.hide)

But I want to pull in and pull out the description window. It should work like example showing below.

Is there any way to do this using Qt designer?


texts: What to do with polygons that were created from a raster file?

texts:
I need some suggestions.  I'm working in Arcmap 10.3.  I used a NDVI landsat product to find patches of row crops.  I then extracted these cells and turned them into shapefiles.  I want to measure the edge of these polygons but am a bit concerned about how the raster format will affect the measurement.  As you can see in the photo, I feel like the zig zag shape will give an over estimation of edge length.

Do you think tools like aggregate and smooth polygons might help me?


texts: How can I make a "smart" infowindow in CartoDB?

texts:
I embed CartoDB visualizations on a variety of websites with different column widths, and on the narrower ones infowindows can be annoying.

The normal behavior on click is for the map to pan to where the infowindow is fully visible above and right of the marker. But if you have a lot of data and you want people to be able to compare multiple markers or polygons, users start getting bounced all over the place.

What I'd like to achieve is something like this Google maps example:

http://gmaps-samples-v3.googlecode.com/svn/trunk/smartinfowindow/smartinfowindow.html

If you click on one of the markers, pretending it's the tip of a clock hand and dragging it clockwise around the frame, you'll see what they're calling a "smart infowindow" -- it repositions itself toward the center of the frame.

But I'm not as interested in the repositioning capability (although that's awesome, it seems like overkill) as I am in it initially detecting which direction to open. If the marker is on the right side, the infowindow pops open to the left. If it's on the left side, it pops to the right.

How can I do that with a CartoDB map?


texts:
I am trying to automate a process of a shapefile. We are using Google Maps engine and we have limitation on number vertices in a shapefile. So I am using a grid to split the shapefiles.

In GUI , I add a grid from research tool then vector->geoprocessing->intersect and get the result but when I try it using API ( with below code ) result is not same. it gives me error as it cannot insert LINESTRING in POLYGON shapefile. this is my code:

QgsApplication.setPrefixPath('C:\\OSGeo4W\\apps\\qgis', True)
QgsApplication.initQgis()
covLayer = QgsVectorLayer('C:\\Shapes\\small coverage.shp','coverage','ogr')
gridLayer = QgsVectorLayer('C:\\Shapes\\small grid.shp','grid','ogr')
qgis.analysis.QgsOverlayAnalyzer().intersection(gridLayer,covLayer,'C:\\Shapes\\intersected.shp',False)
QgsApplication.exitQgis()

Any insight?


texts: PyQGIS QgsVectorLayer.uniqueValues() override active object filtering

texts:
This code snippet gives me the unique values for a vector layers attribute called myfieldname

l = iface.activeLayer()
attrIndex = l.fieldNameIndex('myfieldname')
uniValues = l.uniqueValues(attrIndex)

With an active object filter I only get the unique values of a subset matching the filter criteria.

Is there a way to override this behaviour, in order to get all unique values for an attribute (of the complete layer) even with an active object filter?

(applies also for featureCount() etc.)


texts: Converting units of m day-1 (ERA5) to mm (CRU)

texts: I have downloaded 2 datasets, one from CRU (https://data.ceda.ac.uk/badc/cru/data/cru_ts/cru_ts_3.00/data) of cru_ts_3_00.1901.2006.pre.nc.gz and other is ERA5 data (https://cds.climate.copernicus.eu/cdsapp#!/dataset/ecv-for-climate-change?tab=form) of Precipitation. The CRU data units for precipitation is mm as stated in https://dap.ceda.ac.uk/badc/cru/doc/File_Format_CRUTS3.00.pdf and ERA5 data is in m day-1 as stated in https://cds.climate.copernicus.eu/cdsapp#!/dataset/ecv-for-climate-change?tab=overview. Now I have opened and converted both of the datasets into raster format in R. However I want to make the units same. How can I convert m day-1(ERA5 data) to mm(CRU data)?


texts:
I need to get all features' names, and for that I did a listview to put all features' names on it.

So when I click one feature, not on the map but on the listview, the feature must be selected.

Here it is my code:

for(var f=0;f<select.features.length;f++) {

        if(select.features[f].attributes.ADI == "TALAS") {
        alert(select.features[0].attributes.ADI);

        break;
        }
}

When I click a feature on a map, the select.features.length will be 1.

But I don't need to select the feature on the map, I must have the array of features and then I search attributes on this array and select that feature. That's what I would like to get.

I appreciate for any hint to get this working. I try to write some codes about that:

select.removeAllFeatures();
filter = new OpenLayers.Filter.Comparison({
                type: OpenLayers.Filter.Comparison.EQUAL_TO,
                property: "ADI",
                value: "TALAS"
});

filterStrategy = new OpenLayers.Strategy.Filter({filter: filter});
groups = new OpenLayers.Control.GetFeature({
            strategies: [new OpenLayers.Strategy.Fixed(),filterStrategy],
            protocol: OpenLayers.Protocol.WFS.fromWMSLayer(ilce)
});

select.addFeatures(groups.feature);

I got this error message: Uncaught TypeError: Cannot set property 'layer' of undefined


texts: I have a test map with some choropleth layers. I want to display a legend for every layer, but I don't know how. I want to display only a simple picture as a layer, with turn on layer, display image.


texts: Creating dummy variables between layers based on their geometries in QGIS

texts:
I have the following data:

administrative boundaries for the United States: https://biogeo.ucdavis.edu/data/gadm3.6/shp/gadm36_USA_shp.zip
vector layer containing coastlineshttps://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_coastline.zip
river layerhttps://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_rivers_lake_centerlines_scale_rank.zip

Question I. For a boundary layer
How to create a dummy variable being 1 if a district lies on the coast and 0 otherwise. Alternatively, I also want to know the distance from the district to the closest coast.

Question II. For a river layer
How to create a dummy variable where the dummy is 1 if a county is crossed by a river and 0 otherwise. The thing now is, that each river is assigned a scale rank variable. I only want the dummy to be 1 if rivers with scale ranks between 1 and 5 cross the county
.

I still did not manage doing the task using the field calculator.


texts: How to get SSURGO Soil Shapefile Directly inside QGIS?

texts: SSURGO Soil data can be downloaded from websoilsurvey https://websoilsurvey.sc.egov.usda.gov/App/HomePage.htm by creating an AOI and then checking out, but this is a tedious process. Is there a way to get the Soil Layer Directly inside QGIS by just giving an area of interest?


texts:
Is there a clever way to visualise more than one property in a single map overlay?

Ideally I would like to produce a dot grid map in which the size of the dot represents a scaled property and the colour range would represent a second property of the same feature.

As an alternative this could be achieved by using clustering.

There is a CartoDB clustering example. However, this example doesn't really show how to apply a scaled size based on an aggregated property. There is also no colour based on an aggregated property.

Maybe someone can nudge me in the right direction to show an aggregated property in a hexagon map. The CartoDB wizard has only an option to display point density.

#admin_soa_pop{
  [points_density <= 0.00000416004431488379] { polygon-fill: #BD0026;  }
  [points_density <= 3.71432528114624e-7] { polygon-fill: #F03B20;  }
  [points_density <= 1.4857301124585e-7] { polygon-fill: #FD8D3C;  }
  [points_density <= 7.42865056229248e-8] { polygon-fill: #FECC5C;  }
  [points_density <= 7.42865056229248e-8] { polygon-fill: #FFFFB2;  }
}


texts: Osm2pgsql - Lines instead of polygons that represent buildings

texts:
So I am using a PostgreSQL database with the PostGIS extension to extract building data with addresses from osm data.

I just discovered that in my dataset, imported with osm2pgsql of the german city hamburg there are some entries in planet_osm_line which have addresses. Those entries are also valid polygons, i tried to transform those LINESTRINGs with ST_MakePolygon and that worked for every entry.

Here is an example of what I am talking about:

So I wondered why osm2pgsql did not interpret them as polygons. To my understanding a valid polygon (closed way) needs identical start and end nodes, which is so far true for every entry.

I might be missing something here, so can anyone with more experience can explain this to me?

Here is the text representation of the geometry from above:


texts: Automatically make maps from point file?

texts:
Please forgive me if this is an overly basic or broad question. I am a competent user of ArcGIS.

I have a shapefile containing thousands of points. I'd like to automatically make a map of each point, in which:

The name of the finished map is the name of an attribute of the point
The point is symbolized a certain way
A polygon is chosen from another shapefile based upon the location of the point
That polygon is symbolized a certain way
Text populates a box on the map containing attributes from both the point and the polygon
The polygon is labelled with an expression and in a certain format
among a few other requirements.

I can batch export mxds to pdfs. I can batch rename documents. I'm familiar with data driven pages. I can make one of the maps in a few minutes, but I have hundreds to make.  Help point me in the right direction - is this a job for Python, Model Builder, or something else? I'm willing to spend some time learning, but I'm not sure where to start.


texts: Re-projecting shapefiles (with different projections) to the projection of a specified shapefile: Error 000628

texts:
I am trying to write a script to make a tool that will run in Arcmap 10.7.1.

I have a file with multiple different shapefiles with multiple different projections. Basically, it has to take the spatial reference of a specified shapefile and reproject it to all the others. When I run the script (hard-coding the 'arcpy.GetParameterAsText()' variables to test the functionality of the script), I get ERROR 000628: Cannot set input into parameter out_coor_system.. How can I get around that and make the script functional?

Here is the code I have:

#Import modules
import arcpy, os

#Set workspace directory
from arcpy import env

#Define workspace
inWorkspace = arcpy.GetParameterAsText(0)
env.workspace = inWorkspace
env.overwriteOutput = True

#Define local shape file to reproject to
targetFeature = arcpy.GetParameterAsText(1)

#Describe the input feature class
inFc = arcpy.Describe(targetFeature)

#Define the spatial reference of the input feature class
sRef = inFc.spatialReference.name

#List the feature classes in the workspace
fcList = arcpy.ListFeatureClasses()

#Reproject all the feature classes:
for fc in fcList:
    desc = arcpy.Describe(fc)
    if desc.spatialReference.name != sRef:
        print "Projection of " + str(fc) + " is " + desc.spatialReference.name + ", so re-defining projection now:"
        newFc = arcpy.Project_management(fc, desc.name + "_projected", sRef)
        arcpy.AddMessage(arcpy.GetMessages())
        newFc = arcpy.Describe(newFc)
        print str(fc) + " has been reprojected to " + str(newFc.spatialReference.name)

I'm trying to add the "_projected" to the original file name, to signify that it has been reprojected. Also, I'm using IDLE as my IDE.


texts:
The editor pane in the python console has a button "Open in External Editor".
How do I tell QGIS what program to use as the external editor?


texts: I can't get soilgrid image with Python - OWSLib

texts:
I'm trying get an imagen GTIFF with OWSLib in Python but I get error. Would you help me?
Code:
from owslib.wcs import WebCoverageService
wcs = WebCoverageService('http://maps.isric.org/mapserv?map=/map/phh2o.map', version='1.0.0')

bbox = (-54.7968690628577, -28.35600162342199, -54.77831702365789, -28.33600341642725)

response = wcs.getCoverage(
    identifier='phh2o_0-5cm_mean',
    crs='urn:ogc:def:crs:EPSG::4326',
    bbox=bbox,
    resx=250, resy=250,
    format='GEOTIFF_INT16')

response.read()

Output:

b'\n<ServiceExceptionReport version="1.2.0"\nxmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wcs/1.0.0/OGC-exception.xsd">\n  msImageCreate(): Image handling error. Attempt to allocate raw image failed, out of memory.\n  \n\n'


texts: Does FME remember processes that have run?

texts:
Just getting into FME after being used to ArcMAP Modelbuilder.

Every process in the model in FME is run everytime you press run. In ArcGIS Model builder it remembers the processes run and only runs the processes required to get to the step you are in that have not been run. This only resets if the inputs change or if you reset the model inputs.

Can FME do the same? I want to use it in research so I am building my process as I go on a large dataset so I don't want to waste time to re-run processes that I know work correctly.


texts: Extracting text within parenthesis in a QGIS attribute table

texts:
I have a shapefile of data for Japan, and I need to extract the English character names in my "name" field in order to label the layer. Much of the data is all in Japanese, but some of the names have an alternate English name which can be found inside a set of parenthesis. I'd like to extract this text to another field and use it for labeling. I'm using QGIS.

Example:

I need to turn

野川サイクリング道路 (Nogawa Cycling Road)

Into

Nogawa Cycling Road

I'm attempting to use the regexp_substr function in the Field Calculator:

regexp_substr( "name", '\((.*?)\)')

When I try to run this, I get an error:

  An error occured while evaluating the calculation string: Invalid regular expression '?(.*?)?': bad repetition syntax.

Oddly, in testing, I tried this same regex in PostGIS and it worked great. Does anyone know of a regex that will do this task in the QGIS field calculator?