texts: How to install GeoServer on Ubuntu 11.10 64 bit?

texts:
I borrowed the user name from friend.

I need to install GeoServer on Ubuntu 11.10 64 bit.

first, I installed Oracle Java 7 (JRE/JDK) successfully.

then, i thought i need to install JAI and ImageIO (jre and jdk amd 64-bit) from here
http://download.java.net/media/jai/builds/release/1_1_3/

http://download.java.net/media/jai-imageio/builds/release/1.1/

following these instructions :
http://docs.geoserver.org/latest/en/user/production/java.html#installing-native-jai-on-linux

but that didn't work.

then tried these instructions (for JAI files) :
http://java.sun.com/products/java-media/jai/INSTALL-1_1_2.html#Installation

install done but with these errors :

Cannot create lib/amd64/libmlib_jai.so

Cannot create lib/ext/jai_core.jar

Cannot create lib/ext/jai_codec.jar

Cannot create lib/ext/mlibwrapper

also when i tried to install IimageIO files i got: Please do not attempt to install this archive file

I am not sure what to do now.

maybe it is better if i use 32 bit machine instead of 64 bit ?
or use openJDK instead of oracle java ?
I really need your suggestions ppl.

thank you


texts: Saving parameters for custom python form?

texts:
I'm using QGIS 2.14.0. I have a custom python script that shows a form and then creates some layers from supplied parameters, but I want to save them on running/closing the form, but cannot figure it out (not default parameters, but persistent data from the last dialog).

I saw another kind of similar question and it said in settings/options tick 'reuse last entered attribute values', but this does nothing. I've a feeling it only works for built in forms.

this is the start of my code, can anyone offer a solution to this?

##NameArea=string
##ShowBoundary=boolean
##ShowRoads=boolean
##ShowStructures=boolean

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
import processing
import os

If the solution is to somehow read/write data to a form I can do this but don't know where to start. I just figured this feature was pretty much a pre-requisite for most forms so has a simple solution.


texts: PyQGIS Composer Map size/bounding box problem SVG

texts:
I expect the following code to find the bounding box around feature (which is the magenta polygon in the picture) and render that so that the height is 1000px. The white background is indeed the bounding box of the feature but it's smaller than 1000px, on a bigger gray background, and the whole SVG is still larger; moreover, features outside the bounding box are drawn.

How do I draw the exact extent I want, at an exact specified extent?

from PyQt4.QtSvg import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *

layers = QgsMapLayerRegistry.instance().mapLayers()
for layer_name, layer in layers.items():
    if 'regions' in layer_name:
        break

feature = next(layer.getFeatures(QgsFeatureRequest(4)))
region_name = feature.attribute('regionname')
region = feature.geometry()
bb = region.boundingBox()

renderer = iface.mapCanvas().mapRenderer()

c = QgsComposition(renderer)
c.setPlotStyle(QgsComposition.Print)

composerMap = QgsComposerMap(c, 0, 0, 1000, 1000)
composerMap.setNewExtent(bb)
c.addItem(composerMap)

print (composerMap.rect().width(),composerMap.rect().height())
print composerMap.extent().toRectF()

c.setPaperSize(composerMap.rect().width(),composerMap.rect().height())

svg = QSvgGenerator()
svg.setFileName("/tmp/mwe.svg")
svg.setSize(QSize(composerMap.rect().width(),composerMap.rect().height()))
svg.setTitle(region_name)

painter = QPainter(svg)
c.render(painter)
painter.end()


texts: I would like to configure a Geoevent Extension service to send an email to a group of users when feature 1 intersects feature 2.  Do these features need to be feature services or do map services suffice?


texts: ArcMap Create Lines from Points based on (From-To Relation)

texts:
in ArcMap, I have (Points.shp) which has around 2000 points as shown:

And I have CSV file (let's call it Relation) which contains FromNode, ToNode (Relation of connecting points from-to)

Btw, The numbers in FromNode & ToNode refers to the id of the Points.shp

What I need is how to create (lines) which connecting these points based on the relation (From-To), using any tool, script, pyscript.


texts: Loading GeoJSON into map using JavaScript API in ESRI

texts:
I'm new to ESRI Map.  I'm working on latest version of JavaScript API.

I'm trying to load GeoJSON from this link into my Map.  I'm able to load GeoJSON from ESRI Feature Service URL.

I wanted to load GeoJSON from an object.  My code goes like this.

require([
      "esri/Map",
      "esri/views/MapView",
      "esri/layers/FeatureLayer",
      "esri/widgets/Legend",
      "dojo/domReady!"
    ], function(
      Map, MapView, FeatureLayer, Legend
    ) {

      var defaultSym = {
        type: "simple-fill", // autocasts as new SimpleFillSymbol()
        outline: { // autocasts as new SimpleLineSymbol()
          color: "lightgray",
          width: 0.5
        }
      };

      var renderer = {
        type: "simple", // autocasts as new SimpleRenderer()
        symbol: defaultSym
      };

      var map = new Map({
        basemap: "streets"
      });

      var view = new MapView({
        container: "viewDiv",
        map: map,
        center: [-85.050200, 33.125524],
        zoom: 4
      });

      var usOutline = ...... A single feature from the feature collection in link.
      var featureLayer = new FeatureLayer({
        source: usOutline.geometry,
        title: 'Hello',
        geometryType: "Polygon",
        renderer: renderer
      });

      map.add(featureLayer);

    });

How can I load the GeoJSON?


texts: In ArcGIS 10, is there a way to add a field to a feature class even if it currently contains data?


texts: How to clip rasterdata for later processing?

texts:
My need to dwell with geographical tools is only temporary, thus I am not well versed and my question is difficult to pinpoint exactly. But here it goes:

-I have about 5-10 sites, for which I need to assemble geographic info from many sources (satellite data, aerial images, meteorological data, geological info, many of these sources having several layers). These come with many formats (JPEG2000, Geotiff,  HDF5, ...)

-I can use either Matlab, python (arcpy?), QGIS or ArcGIS.

-The essential part is to drop the extra data from various files, to clip the raster to contain only the small rectangular area (the location is known) and write this down (likely in ASCII format for debugging purposes). How to do this with some of the tools? The process is not repetitive. (This is just for practical value: instead of 2G of data I will have only 20Mb or so.)

-the next phase establishes a new raster grid and interpolates the data to those grid points. It is nice if this can be done with the tools but quite likely I have to do it
in the script, since the grid point data vectors are later fed into other data analysis software.

The grid sizes vary from 2km to 2 meters. Is one of the tools enough?
Should I use several tools for addressing different formats?
Is there reason to learn the tools or just use one week to code a script in Python?


texts: Changing color rendering of multiple raster layers at once

texts:
I have many raster layers in my QGIS project and they're all currently rendered as:
Multiband color
Red band: Band 1
Green band: Band 2
Blue band: Band 3
I'd like to know if there's any way I can change the rendering of all raster layers to:
Red band: Band 3
Green band: Band 2
Blue band: Band 1
without having to do it by hand, layer by layer.


texts: I would like to know if it is possible to use your own map created with Google Maps Engine with OpenLayers? Any examples?


texts: Interpolate biweekly NDVI to daily using Google earth engine

texts:
How to interpolate time series data using GEE?

For example, the Landsat have 2-3 observations per month, I want to interpolate the series of NDVI to daily scale and get 365 observations a year, then I can use daily series data to calculate vegetation phenology information e.g., the onset of growing season.

I want to know if there is some function to realize this.


texts: Display all features of layer using GetFeatureInfo in leaflet

texts:
I wan't to display all the features of the WMS layer using GetFeatureInfo in leaflet.The features like point,state name or label and when I click on point I wan't to open new HTML page and display features without popup from Geo-Server.
I have done this code but output showing blank map window.Thanks

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

About

Urban Expansion

Cities

Data

Historical Data

Contact


texts: Why is the QGIS map composer grid showing unreasonable lat/lon?

texts: I need to print a map with coordinates in the format of Degree, Minute and second. CRS is set to WGS 84 (EPSG 3857) for both layers and on the project properties it is set to enable transformation on the fly. But then, when in map composer I get the following coordinates, which make no sense to me.


texts: Renaming multiple rasters with ModelBuilder gives errors 000622 and 000878?

texts:
I have 101 raster images stacked inside a folder. The Rasters are of arcinfo grid type (*.asc) while having  a long string name 'NW2D_T100_05_base_depthXXX.asc' where XXX ranges from  000 to 100. What I would like to do is to add each of them with a single DTM. I tried to use ArcGIS 10.1 model builder by inserting raster iterators and raster calculator tool and the model looks as follows:

When I ran the model, the following error message was displayed:

The error numbers are 000622 and 000878.

It looks like that I have a very long string name for each raster data. How can I rename the raster with the 8 digits of the name (i.e., depthXXX.asc) for each file ? And how can I save 101 new rasters in  anew directory.


texts:
My problem is, I tried this solution in order to clear the polylines on the map. Here's some code for your reference:

var line;
 $scope.clearPolylines = function(){
    for(var i = 0; i < coordinates.length; i ++){ // coordinates is an array
        map.removeLayer(line._latlngs[i]);
        console.log("[ DELETE LINES ]", line._latlngs[i]);
    }
}

But it is not clearing on the map. Are there other ways possible?

This solution is based on  https://stackoverflow.com/questions/14585688/clear-all-polylines-from-leaflet-map


texts:
I am trying to add a GeoJSON file to my OpenLayers application and I am getting a cross origin request error.

My main question is can a GeoJSON be stored locally (within a folder on my C drive) and be used in OpenLayers?

I am just learning OL and I just want to be able to test some of the capabilities without setting up a web server.

Secondly, if its not possible to bring in a GeoJSON without having a web server or something like jetty running on my computer, is there a better format for bringing in features to OpenLayers?

Are KML files a better option for what I am trying to do?


texts: Google Maps click event does not work on mobile devices

texts:
I created a map using Google Maps API with Angular Material libraries (https://material.angularjs.org/)

After I upgraded Angular Material from 0.70 to 0.82, I found the info window did not show up by clicking the map on mobile phones or tablets.It looks like the map click event cannot be fired on touch screens. This issue did not happen when using the version 0.70.

The old map w/ Angular Material 0.70:
http://csj-councildistricts.appspot.com/index.html

The new map w/ Angular Material 0.82:
http://6.csj-councildistricts.appspot.com/index.html

Does anyone have the similar issue?


texts: Querying multiple shapefiles using Python and GDAL:OGR

texts:
I would like to select everything from one shapefile based on matching ID's in another shapefile. But I am unsure about how to proceed with the syntax using Python.

Mainly because .ExecuteSQL('sqlstring') is applied to a single data source, after that is loaded with the driver.Open('filePath') command. So something like this below doesn't function:

inFile0 = r'C:\data\buildings\buildingPoints.shp'
inData0 = driver.Open(inFile0, 1)
inLayer0 = inData0.GetLayer()

inFile1 = r'C:\data\buildings\buildingPolygons.shp'
inData1 = driver.Open(inFile1, 1)
inLayer1 = inData1.GetLayer()

sql = "SELECT * FROM buildingPolygons WHERE FEAT_ID IN (SELECT FEAT_ID FROM buildingPoints WHERE NET = %s)" % sampleID

samplePolygons = inData1.ExecuteSQL(sql)

When I create two separate sql statements inside a script running sequentially, where the first selects the relevant points from the points dataset and the second selects the building polygons based on matching the IDs with the selected building points, the result is empty:

sql0 = "SELECT" * FROM buildingPoints WHERE NET = %s" % sampleID
samplePoints = inData0.ExecuteSQL(sql0)

sql1 = "SELECT" buildingPolygons.* FROM buildingPolygons, samplePoints, WHERE buildingPolygons.FEAT_ID = samplePoints.FEAT_ID"
samplePolygons = inData1.ExecuteSQL(sql1)

How could I make this work? Or should I look into using the subprocess module and run ogr2ogr using as input a .vrt of the two shapefiles and apply the sql there?


texts:
When I try create a new database, Spatialite plugin return the dialog error box:

The SQL query seems to be invalid, no such table: views_goemetry_columns

When I click in OK button, a second dialog box (Transform database into SpatiaLite) is open with the text :

This database seems to be a valid SQLite database but not a valid SpatiaLite One
Would you like QSpatiaLite to automatically transform it to SpatiaLite (information won't be lost) ?

When I click in the YES button, a new dialog error box open with the text:

Unable to convert Database to SpatiaLite

When I click in OK button, another dialog error box open with the text:

The SQL query seems to be invalid.

The database was set up, but I can't import spatial layers.
How resolve this problem?


texts: ToolStrip Button click event is not fired on first click

texts: I have a dockable window in a custom ArcMap extension.  The dockable window has a regular windows forms ToolStrip control in it with ToolStripButtons on the ToolStrip.  When the parent dockable window is not Active and I try to click a button on the toolstrip within the dockable window, I have to click it twice - once to set the focus and then again to fire the click event.  I only want to click the button once even if the dockable window is not currently active.


texts:
I am trying to do a research project on changes in Land Cover in a London borough over time. I have access to OS MasterMap data through Digimap (the service that provides this data to universities in the UK), but I can't see a way to download a time series of data directly.

I've noticed that the tables that are produced from the OS MasterMap data have a couple of fields that give the date the feature was changed, and whether the feature was New or Updated at that point. However, I can't see a way to use this data to produce previous versions of the map data. Is this possible to do? It would seem strange for them to include these fields if there was no way to process the data to make them useful, but I can't see anywhere that they say (for example) what was there beforehand, if this feature is new.

I have heard that the British Library have snapshots of OS MasterMap data from various times, but these are only available as images, not as GIS data.

Does anyone have any ideas as to how I could obtain a time series?


texts: Problems with calling GeoJson properties in OpenLayers3 (v3.1.1.)

texts:
I have a problem, property from GeoJson is not recognizable by OpenLayers3.
Debbuger says: Uncaught ReferenceError: feature is not defined on line futher below in OpenLayers3 code:

var icon = feature.get("icon");

GeoJson

{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "geometry": {
                "type": "Point",
                "coordinates": [
                    12.231,
                    30.342
                ]
            },
            "properties": {
                "temperature": "21.1",
                "icon": "4"
            }
        },...

OpenLayers3 v.3.1.1.

var icon = feature.get("icon");
var styleTemperature = new ol.style.Style({
        image: new ol.style.Icon({
        src: 'http://prognoza.hr/metsimboli/' + icon + '.gif'
})
      });

var weatherLayer=new ol.layer.Vector({
        layer:'weather',
        visible: false,
        source: new ol.source.GeoJSON ({
                   projection: 'EPSG:4326',
                   url: 'path/to/geojson.geojson'
        }),
        style: styleTemperature
});

Can somebody suggest me something?

EDIT:

This part of the code works, but there is no added style to it, or should I say style is default:

var weatherLayer = new ol.layer.Vector({
        layer:'weather',
        source: new ol.source.GeoJSON ({
                projection: 'EPSG:3765',
                url: 'path/to/file.geojson'
        })
});


texts: Arcgis Pro's Eliminate not eliminating selected polygons?

texts:
I'm trying to use Eliminate to reduce the number of polygons in my dataset that have a very small area size. As far as i understand the tool, you have to first create a selection in the dataset, and afterwards run the Eliminate tool. What didn't make sense for me, was the fact that the number of selected features didnt add up with the number of features that were returned in the elimination process. I also tried to calculate the area for each of the newly created features, and there were still some as small as the selection i made previously.

below are shown the first selection (47676) and the expected number of returned rows is shown (18.527 vs 20395).
I tried to run the elimination a second time, and the returned rows were still not adding up (18531 vs 18640).

I read the part that selected polygons that are touching other selected polygons, wont be eliminated. This is however not always the reason for the tool not eliminating parts.

IN this example, the yellow polygon is fully contained by the turquise/grey polygon, and should have been eliminated by the tool. All polygons with a value in POLY_AREA smaller than 100 was supposed to be eliminated.

I was therefore wondering if anyknody can explain what is happening here?


texts: Calculating weighted values of field in Attribute table in QGIS

texts:
In QGIS in the Attribute Table, I have a field with values.

For this field, I need to calculate the Sum of all its values and divide by the number of values that are not equal to "NULL".

Afterwards, I need to create a new field, where each value of the original field will be divided into the previously calculated (the Sum of all its values and divided by the number of values).

Is it possible to do this with the function editor in the Field Calculator (as this action will be repeated more than once)?


texts: How can I do feature selection on WMS layer (added JSFiddle example)

texts:
I have WMS layer with roads serverd from Geoserver (MS SQL View Layer)

http://jsfiddle.net/Alophind/5pSrJ/

I want the user to be able to click on a line and "highlight" it , set it selected so he can use edit tools (split,union) to edit it.

How can it be done ?


texts: Seeking software for opening .img files?

texts:
I have some cpc precipitation data stored in IMG format.

Could anyone suggest a software beside ERDAS imagine or Arcgis which can be freely downloaded?


texts: How to set default bad layers handler in the python plugin?

texts:
I am trying to write a python plugin for QGis. My python plugin will copy a template qgis project into data catalog with ESRI shapefiles. Because data catalogs can vary in its content there will be possibility that in the output catalog may be less layers than in the project template.
Is it possible to set default bad layer handler inside python plugin every time when I load programmatically copied project form output catalog?
I am new in programming qgis plugin and I suppose that I should use something like that:

proj = QgsProject.instance()
handler = proj.setBadLayerHandler() #I don't know what kind of arguments I should put in this method.

Thanks,


texts:
Relating to a previous post, found an excellent source for data called Koordinates.

I downloaded this 5m DEM of New Zealand, Wellington.
This is in a TIFF format and I have added to ArcMap 10, applied a colour ramp, some hillshading, and applied a standard deviation on the histogram:

I know that the MSD format does not support hillshading, so I published direct from the MXD.  As a dynamic service, it draws fine (but slow as you would expect).

So I cached the service using Mixed mode, and the compact storage type.

Viewing the cache in both ArcGIS Desktop and REST, the resulting cache is a solid red:

At some of the closer levels, the expected symbolgy comes through in some tiles:

Have tried re-creating several times and still getting the same problem.

Would like some advice on troubleshooting this.
Have not seen anything related in the Esri forums.


texts: Spatial databases in location based search

texts:
I'm developing a location based application in which i was in the plan of dividing the globe into grids and when the client requests information, will create a table or similar data model which will expire in 5 mins and serving the same data for the clients within the same grid within that period. Which i explained in the following StackOverflow question.

https://stackoverflow.com/questions/19667637/caching-location-based-data-cluster

But i come across the term spatial databases today and now i like to know how it differs from the traditional model when searching from the databases of million records especially in terms of speed? Will it really speed up data retrieval of location based data? If so, can you explain how?

Problem: The problem is, i was in the idea of using traditional database, but if i use traditional DB and if i have more than a million records, when we execute each search query it needs to search through the the database of million records, so i planned to,

Divide the globe by some sort of grid each of around 1 sq.km and when we showed all the values from that grid,then the next nearby grid and so on which i need to save in some buffer.
And when someone request the data from the same grid in the next 5 mins, i don't want to search the entire DB again, and i need to give the values from the buffer and not querying the entire DB again.

Here the main problem is, buffering the data and retrieval from it to speed up the process and minimizing the server load. So i like to know that this spatial DB simplify these process in any way? Or is there any DB already available which provides any grid system like this?

@ike: What i meant as speed is the performance, i.e., the number of rows returned against the number of rows on which the query runs on(may be the entire table of few million records). I'm trying to find some way to reduce the number of rows(or the entire table) on which the query runs


texts: QGIS Problem in calculating the distance of points from line

texts:
I need to calculate the distance of villages ( villages as points) from a border. I have a shapefile containing the locations of villages and another one containing the border. On QGIS, I ran the 'distance to nearest hub' command to calculate the distances. The distances that the command returned looks absurd. Can anyone help to locate where potentially am I going wrong?
(I have checked that both the layers are in the same coordinate reference system.)


texts: ArcGIS Server: Limit creation of cache tiles to certain areas by providing a feature class

texts:
I am publishing a map service to ArcGIS Server 10.0. I cannot succeed to limit the creation of tiles to certain areas only, instead of creating them for the entire map. According to this help document ("Map caching based on feature boundaries
"), as far as I understand, this should be possible.

In my map document I have an aerial photo of an area containing some villages. I have a shape file polygon feature class with rectangular polygons containing the villages. I want to provide this feature class to limit the creation of tiles only to the portions of the aerial photo within these bounding rectangles. I am not interested in the rural areas in between the villages. It would be a waste of computing time and disk space to have tiles created outside the rectangles.

I followed the steps in the above mentioned help document closely, providing my shape file here:
Tool "Manage Map Server Cache Tiles" -> Update specific areas using a feature class ->
Update Extents Feature Class (optional) -> MyBoundingBoxes.shp

Unfortunately my feature class is completely ignored. When I examin the resulting PNG tiles, I notice that tiles were created for the entire map, also for the rural areas in between.

Have I done something wrong?


texts:
I have a point layer and need to draw a line passing through or in the vicinity of the points (up to 0.5m ditance between the point and the line).
The nodes of the line will recieve the attributes of the points closest to them.


texts: I have a raster with different values (negative and positive).  I want to color the raster based on a rule:  All positive values have a certain color and all negative another.  I don't want to use 'singleband pseudocolor' becuase it creates several values in the legend.  I just need two values to show up in the legend (red for negative, blue for positive).


texts:
Just set up my MapServer on Ubuntu - all good. Downloaded this tutorial and it works fine: http://download.osgeo.org/mapserver/docs/mapserver-tutorial.zip

I've tried adding a PostGIS layer to the map, without success - can anyone help me. Here is the map file:

MAP
  IMAGETYPE      PNG
  CONFIG "MS_ERRORFILE" "/tmp/ms_error.txt"
  DEBUG 5
  EXTENT         11.8496964862307 55.5990839999999 13.02463 57.8704279999999
  SIZE           400 300
  SHAPEPATH      "../data"
  IMAGECOLOR     255 255 255

  # Layer objects are defined beneath the map object.  You need at least one
  # layer defined in your map file before you can display a map...  You can
  # define as many layers as you'd like although a limit is typically hard-coded
  # in map.h in the MapServer source.  The default limit is set at 100.  You'd
  # have to have a very specialized application to need more than 100 layers in
  # your application.

  # Start of LAYER DEFINITIONS ---------------------------------------------
  LAYER # States polygon layer begins here
    NAME         states
    DATA         states_ugl
    STATUS       OFF
    TYPE         POLYGON

    # The class object is defined within the layer object.  You can define as
    # many classes as you need (well, there are limits as with layers, but it's
    # senseless to define more than ten on a "normal" layer.  There are
    # situations, however, where you might have to do it.)
    CLASS
      NAME       "The Upper Great Lakes States"

      # There are styles in a class, just like there are classes in a layer,
      # just like there are layers in a map.  You can define multiple styles in
      # a class just as you can define multiple classes in a layer and multiple
      # layers in a map.
      STYLE
        COLOR        232 232 232
        OUTLINECOLOR 32 32 32
      END
    END
  END # States polygon layer ends here
  LAYER
    CONNECTIONTYPE postgis
    NAME "plats"
    DEBUG 5
    # Connect to a spatial database
    CONNECTION "user=postgres dbname=test host=localhost password=xxxxx"
    #PROCESSING "CLOSE_CONNECTION=DEFER"
    # Get the lines from the 'geom' column of the 'plats' table
    DATA "geom from plats"
    STATUS OFF
    TYPE POLYGON
    CLASS
      NAME       "Plats"

      # There are styles in a class, just like there are classes in a layer,
      # just like there are layers in a map.  You can define multiple styles in
      # a class just as you can define multiple classes in a layer and multiple
      # layers in a map.
      STYLE
        COLOR        232 232 232
        OUTLINECOLOR 32 32 32
      END
    END
  END
  # End of LAYER DEFINITIONS -------------------------------

END # All map files must come to an end just as all other things must come to...

the debug output wasn't much help ...

[Fri Mar 14 19:28:07 2014].610665 CGI Request 1 on process 8022
[Fri Mar 14 19:28:07 2014].611179 msDrawMap(): rendering using outputformat named png (AGG/PNG).
[Fri Mar 14 19:28:07 2014].611200 msDrawMap(): WMS/WFS set-up and query, 0.000s
[Fri Mar 14 19:28:07 2014].611318 msDrawMap(): Layer 0 (states), 0.000s
[Fri Mar 14 19:28:07 2014].611336 msDrawMap(): Drawing Label Cache, 0.000s
[Fri Mar 14 19:28:07 2014].611343 msDrawMap() total time: 0.001s
[Fri Mar 14 19:28:07 2014].616884 msSaveImage(stdout) total time: 0.006s
[Fri Mar 14 19:28:07 2014].616976 mapserv request processing time (msLoadMap not incl.): 0.006s
[Fri Mar 14 19:28:07 2014].616989 msFreeMap(): freeing map at 0x8211eb8.
[Fri Mar 14 19:28:07 2014].617015 freeLayer(): freeing layer at 0x821a130.
[Fri Mar 14 19:28:07 2014].617025 msPostGISLayerIsOpen called.

Can anyone help?


texts: Arcgis SDK for Android - Query multiple feature layers

texts:
I'm developing a android application with arcgis sdk.
I have a basemap and then 4 different feature layers on top of it.
I began with the arcgis sample "attribute editor" and removed the edit functioncs because I just want to view the features of the clicked point.

I had a little issue with querying a point or a 1px line on the touchscreen but I solved it with creating an envelope around the clicked point and then query the envelope.

So far so good! But my problem now is, that I can only query one feature layer at a time. I want to click on the map and then it queryies the layer 1, layer 2, layer 3, and finally layer 4. How can you do this? I was just able to query layer 1 OR 2 OR 3 OR 4...

I have this code sections:

featureLayerDefault.selectFeatures(query,
ArcGISFeatureLayer.SELECTION_METHOD.NEW,
new CallbackListener<FeatureSet>() {
    // handle any errors
    public void onError(Throwable e) {
    Log.d(TAG, "Select Features Error" + e.getLocalizedMessage());
    }

if (queryResults.getGraphics().length > 0) {
Log.d(TAG, "Feature found id=" + queryResults.getGraphics()[0].getAttributeValue(featureLayerDefault.getObjectIdField()));
// set new data and notify adapter that data has changed
listAdapter.setFeatureSet(queryResults);
listAdapter.notifyDataSetChanged();


texts: Is PPDM Lite still available as Geodatabase template?

texts:
I have been looking for a Geodatabase template that conforms to the PPDM Lite standard and believe that one was previously available at PPDM.org.  However, that page no longer exists and no Geodatabase download is mentioned at the PPDM Lite 1.1 page.
Does anyone know if one is still available anywhere?


texts:
The arcpy snippet below is intended to generate a list of technicians scheduled to work, and then generates a unique feature layer for each one's stops.

fcStops = "TaskRoute1016"
with arcpy.da.SearchCursor(fcStops, "AssignedEngineers") as cursor:
    TechList = sorted({row[0] for row in cursor})

print "Technicians today:"
for Technician in TechList:
    subsetSQL = "AssignedEngineers = '{}'".format(Technician)
    TechLayer = "Tech {}".format(Technician)
    arcpy.MakeFeatureLayer_management(fcStops, TechLayer, subsetSQL)
    print("  {}: {} stops".format(Technician, arcpy.GetCount_management(TechLayer).getOutput(0)))

It works perfectly once (selects just his stops, makes a feature layer, prints out the count), but only once. After the first iteration, I get:

Runtime error  Traceback (most recent call last):
File "<string>", line 5, in <module>
File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6043, in MakeFeatureLayer
 raise e ExecuteError:
ERROR 000622: Failed to execute (Make Feature Layer). Parameters are not valid.
ERROR 000628: Cannot set input into parameter in_features.

The 000628 error says "The syntax used for the parameter is invalid for the parameter type" -- I'm not sure what's changing in my syntax, however, since the in_features variable should still be fcStops each time through the loop. If it's working once, why not multiple times?

I am running this in the Python window in ArcMap 10.2.1.


texts: How to crop TIFF image without losing classes in QGIS 3.2?

texts:
I am trying to crop a TIFF file in QGIS 3.2. What I am currently doing is going to Raster -> Extraction -> Clip Raster by Extent... . Then I type in my min and max x,y values. The new raster layer is the cropped image but all of my layers from the original TIFF image are gone. Attached is a screenshot of the before (noncrop.jpg) and after (crop.jpg) images. Note in the photos the layers are different for the Cresults raster vs. the Clipped (extent). How do I clip the Cresults TIFF image without losing its classes?

Figure 1. QGIS screenshot of the original TIFF image with multiple classes.

Figure 2. QGIS screenshot of the original TIFF image with the Cropped image overlayed. Note the missing classes in the cropped image.


texts: Display data as georeferenced using three matrices of lat / long / data?

texts:
I have a 12 x 12 matrix of data and two more matrices of corresponding latitude and longitude.

Is there a way I can import the data as a georeferenced raster in ArcMap?

For example, in Matlab, plotting pcolor(lon,lat,data) would project the data correctly. I have tried and failed to save the file as a geotiff in other software, so I must try now to display it in this way.


texts: Getting my data to be focused on one state

texts:
I have a map I'm working on for a final project and I'm a bit stuck as what to do moving forward.

I have a map I'm making for forest fires and the type of tree in that area and I'm trying to get it zoned in on Colorado. I have the shapefile for Colorado in my map and I have the forest fire information, and I even have a map of the US for the type of tree in each location.

But what I'm stuck at is trying to get my map to only be focused on Colorado so when I calculate the type of tree on each data point it will only be for Colorado.

Sorry if this is a bit convoluted but I'm just a bit confused and stuck right now. The tree information on my map is an .img file. I'm using ArcGIS 10.7.1.


texts: "undo" pivot table with ArcGIS ModelBuilder?

texts:
I have a table in a geodatabase with an ID and a text field for each day of a week. The text fields contains 0 or 1 for monday, 0 or 2 for tuesday, and so on. One ID can have multiple day entries, so for every day, or only on weekends,...

I need to join this information to another table. It contains a Name, and the ID of the first table. At the end, it should have the following fields: Name, ID, Day. and the day information are coming from the seven "day-fields" and if one ID has multiple day entries it has to be duplicated.

Is there a possibility to do this, or do I need to create seven distinct tables and re-assemble them afterwards?

I am using ArcGIS 10.1 for Desktop.


texts: QGIS Geopackage export of layer symbology

texts:
I want to save several layer as one geopackage. That works fine. My problem is, that the layers in the geopackage have a different symbology than they had before. I found a video where you could easily keep the smybology, but in my QGIS version (3.12.3 Bucaresti) it is not possible.
Is there another way to export a layer as a geopackage file and keep the same symbology?


texts:
I am new to QGIS, it there a way to generate STereoLithography (STL) files with custom output size (I want to match DEM's layer size).

Here it includes 'white' canvas inside model.

I can't change size, it won't let me.


texts:
At my work we use ArcPad to collect inspection point data. One of the fields is SITEID that is manually entered by a staff member (its a combination of their initials and a number (eg ST01). I would like to add an event handler (I'm using OnChange in my code) to my APL file that checks if the SITEID already exists. Im finding it very difficult to debug in ArcStudio so any help would be greatly appreciated!

Below is my APL file

<?xml version="1.0" encoding="UTF-8"?>
    <ArcPad>
<LAYER name="SoilSurvey_InspectionSitesSimple_UTM10N" transparency="1">
    <FORMS>
        <EDITFORM name="EDITFORM" caption="Soil Inspection Sites" width="160" height="150" color="Blue" picturepagevisible="true" attributespagevisible="true" symbologypagevisible="true" geographypagevisible="true" required="false">
            <PAGE name="SiteDetails" caption="Site Details">
                <LABEL name="slblSiteID" x="3" y="3" width="35" height="10" caption="Site ID" tooltip="" group="true" border="false"/>
                <LABEL name="slblDate" x="3" y="18" width="35" height="10" caption="Date" tooltip="" group="true" border="false"/>
                <EDIT name="txtSiteID" x="45" y="1" width="75" height="12" defaultvalue="" **onchange="Call CheckSITEID"** tooltip="" tabstop="true" border="true" required="true" field="SITEID" uppercase="true" multiline="true"/>
                <DATETIME name="dtDate" x="45" y="17" width="75" height="12" defaultvalue="Now ()" tooltip="" tabstop="true" border="true" field="DATE_"/>
                <LABEL name="slblComments" x="3" y="35" width="35" height="10" caption="Comments" tooltip="" group="true" border="false"/>
                <EDIT name="txtComments" x="3" y="48" width="100" height="50" defaultvalue="" tooltip="" tabstop="true" border="true" field="Comments" uppercase="true" multiline="true"/>
            </PAGE>
        </EDITFORM>
    </FORMS>
    <SYMBOLOGY>
        <SIMPLELABELRENDERER visible="false" field="" rotationfield="" expression="" language="">
            <TEXTSYMBOL fontcolor="Black" font="Arial" fontsize="8" vertalignment="bottom" rtl="false" fontstyle="regular"/>
        </SIMPLELABELRENDERER>
        <SIMPLERENDERER>
            <SIMPLEMARKERSYMBOL color="76,230,0" width="3" outlinewidth="1"/>
        </SIMPLERENDERER>
    </SYMBOLOGY>
    <METADATA/>
    <QUERY where=""/>
    <FIELDHISTORY>
        <FIELDS>
            <FIELD name="X"/>
            <FIELD name="Y"/>
            <FIELD name="SITEID" value="3"/>
            <FIELD name="Z"/>
            <FIELD name="DATE_"/>
            <FIELD name="ID"/>
            <FIELD name="Comments" value=""/>
        </FIELDS>
    </FIELDHISTORY>
    </LAYER>

<SCRIPT language="VBScript">
    Sub CheckSITEID
    Dim objEditForm, InputSite, ExistingRec
    Set  InputSite = ThisEvent.Object.Value
    Dim objRS
    Set objRS = Map.Selectionlayer.Records
    If (objRS.RecordCount > 0) Then
        objRS.Movefirst
        While Not objRS.EOF
            Set ExistingRec = objRS.Fields("SITEID").Value
            If InputSite = ExistingRec Then
            ThisEvent.Result = False
            ThisEvent.MessageText  = "SITEID ALREADY EXISTS"
            End If
            objRS.MoveNext
         Wend
    End If
    End Sub
</SCRIPT language="VBScript">
</ArcPad>


texts: How to report that a Sentinel-2 Level-1C granule is not georeferenced accurately (Google Earth Engine)

texts:
I made some greenest pixel composites in GEE and they looked odd--like a river showing up in a town. I found that a granule from 2015 is not georeferenced accurately...then I found another from the same date.

https://code.earthengine.google.com/006fba9fa62c4b9912d38744c9a19d65

Is there a process for reporting this sort of thing?


texts:
I'm building a pretty complex function to rasterise some parameters of my point cloud using grid_metrics (lidR version 2.2.2).

no_fun = function(x,y,z,cl){
  require('lidR')
  require('dplyr')

  metrics = list(
    Prom_class = integer(),
    Tree_ht = numeric(),
    TreeCov = numeric())

  if(length( z[z > 0]) <= 5){return(metrics)}

  [...]
}

(The four variables x,y,z and cl are meant to recreate an internal LAS for an easier handling of the points through the process).

Running the following code:

ctg = readLAScatalog(file.path(wdir,'test.laz')
opt_chunk_size(ctg) <- 30 # to limit the size of an intermediate clustering step
opt_chunk_buffer(ctg) <- 0
opt_select(ctg) <- "xyzc"
opt_filter(ctg) <- "-keep_scan_angle -90 90 -drop_z_above 55"

fuel_metrics = grid_metrics(ctg, ~no_fun(X,Y,Z,Classification), res=10)

I get the aforementioned error:

chunk <- readRDS("\Temp\Rtmpgjmsxy/chunk176.rds")
 las <- readLAS(chunk)
nrow(coords) > 0 is not TRUE

I filled the whole function with a bunch of "warnings" in order to avoid the situation of having less than 5 points and, checking the file, it happens to have almost 20 points (10 of which >0)...

On the web I saw it might be related to the raster package but I don't have any rasterization within the function. I think it might be related to the rasterisation of the outputs and, hence, maybe related with assignment of empty variables (i.e. integer() VS NULL; the former choice seems necessary to avoid another error where the first empty variable's format has to be assigned somehow).

UPDATE 1

Below a full reproducible example with lidR's own data:

no_fun = function(x,y,z,cl){
  require('lidR')
  require('dplyr')

  metrics = list(
    Prom_class = integer(),
    TreeCov = numeric())

  if(length( z[z > 0]) <= 5){return(metrics)}

  # create las
  inplas = data.frame(X=x, Y=y, Z=z, Classification=cl)
  inplas = LAS(data=inplas, proj4string=sp::CRS(projectCRS))
  inplas = lasrescale(inplas, xscale = 0.01, yscale = 0.01, zscale = 0.01) # fix accordingly with the original one

  inplas = lasfilter(inplas, Z > 0)

  if(is.empty(inplas) | nrow(inplas@data)<5) return(metrics)
  crs(inplas) = projectCRS

  las = lasfilter(inplas, Classification != 2)
  if(is.empty(las) | nrow(las@data)<5) return(metrics)

  # initialise output variables
  cl = integer()
  TC = numeric()

  ### other stuff happening here

  # final output
  metrics = list(
    Prom_class = cl,
    TreeCov = TC)

  return(metrics)
}

library(lidR)

LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
ctg = readLAScatalog(LASfile)
opt_chunk_size(ctg) <- 50
opt_select(ctg) <- "xyzc"
opt_filter(ctg) <- "-keep_scan_angle -90 90 -drop_z_above 55"

projectCRS = '+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0'

fuel_metrics = grid_metrics(ctg, ~no_fun(X,Y,Z,Classification), res=10)

UPDATE 2

a workaround is to set empty variables as follows:

  metrics = list(
    Prom_class = as.integer(NA),
    TreeCov = as.numeric(NA))


texts: Open source WMS that supports spatial aggregation and image tiling

texts:
I have daily raster data in geotiff/netcdf format for the world in 10x10 lat-long degrees tiles. Every day, new data is coming in. I want to let people access this data using WMS through a GIS portal.

Ideally, I want to use open source solutions to do so, and therefore would like to make use of e.g. thredds data server. However, Rich Signell pointed out that thredds doesn't support spatial aggregation and image tiling.

So my question is:

What open source WMS software is out there that support file aggregation (multiple files to cover spatially the world) as well as temporal (daily new data)?

Current structure of the data:

I would like to keep the 10x10degrees data tiles, so to cover the world I would have ~648 (360*180)/(10*10) netcdf files. I'm only interested in terrestrial area so realistically there will be less than ~200 netcdf files, each with a grid resolution of 250meter.

I don't want to merge these tiles to create one large file for the entire world, so spatially I have many netcdf files. The reason for this is that one file to cover the entire world becomes too large.

Temporal I would have many files as each day new data comes in for the world. Right now I store this data in separate temporal netcdf files. There is no need to keep this data in temporal separate files for each tile, I could every day write to an existing data netcdf file and add the current data.

In the end I would like to have data presented in a GIS portal such that people can browse to an area of interest and browse through the different time layers using a slider or by entering start and end data.

Is it feasible what I would like to do with anything else than ArcGIS database?


texts: MapGuide with RDBMS database design guidelines

texts:
I'm brand new to GIS programming and I am designing a GIS application and I think I want to use MapGuide because of its user selection features (Digitizing and Redlining - or can I do this easily with MapServer and OpenLayers?). I am reading through documentation/tutorials and I can't seem to find any guidelines to designing a spatial database for use with MapGuide.

Can I design it in any way as long as I keep query optimization in mind? It looks this way from the diagrams, but I want to be sure.

What should I store in my separate DB and what in the MapGuides server repository in .shp files?

Edit: MapGuide stores

Thanks


texts: Create Feature Class does not have any input fields present?

texts:
I'm using ArcMap10.2 and when tying to create feature classes, empty tool is displayed.

when i click OK "error 000735 :feature class:name is required "is displayed.


texts:
I used Mapbasic to create tables. I create only table's structure without any objects inside.
How i want make styles for tables. I tried something like:

Dim b1 As Brush,
p1 As Pen
b1 = Makebrush (1, 0, 16777215)
p1 = MakePen(1, 2, 16776960)

select * from myTable
into __TO__UPDATE NoSelect

Update __TO__UPDATE
Set OBJ = ChangeColor(OBJ, p1, b1)

Function ChangeColor( ByVal oRegion As Object
            , ByVal penNew As Pen
            , ByVal brsNew As Brush) As Object

  Alter Object oRegion
  Info OBJ_INFO_BRUSH, brsNew
  Alter Object oRegion
  Info OBJ_INFO_PEN, penNew

  ChangeColor = oRegion

End Function

But styles still defaults.
After i tried to add point to table and used MakeSymbol() to add symbol style, and its worked.

So MapInfo idialogy means that i gonna have a objects in table to create styles or i just do something wrong?


texts: Find changes between maps using only features

texts:
I have 2 maps with district borders of a region. Each district has its own unique ID number.

From one map to another some districts gave a small part of their territory to the neighboring districts and I want to make a map containing only these small areas that changed districts.

The borders of the 2 maps are badly designed and I can't overlay the maps and pick the areas where the borders change because in the end I would have many areas that come form the bad design and not only the ones that actually change district.

But I know exactly which district lost a part and which district received it, using their ID number.

i.e. An area moved from district with ID:1 to district with ID:2 but I know nothing else about this area.

Is there a way to make a map from these area using only their unique ID?


texts: ArcGIS Pro map units locked as Degrees - how to change to Meters

texts:
Instructions here don't seem to add up. I'm not able to change the map units to Meters.

Is there another way?


texts: How connect GeoServer to existing Mongo database?

texts:
I'm trying to connect existing database in MongoDB to GeoServer.
I connected Mongo to GeoServer, but I can't see exisiting DB. (which I created early by cmd)
However, if I create new db from GeoServer, I can see it in Mongo and in GeoServer.

How published exsisting db in GeoServer?
If I need a special schema, how create it?


texts: How to update a field with a UUID if it doesn't have one already?

texts:
I have a text field called GIS_ID which I would like to populate with a unique UUID only if it doesn't have a valid UUID already. That is, if the field is NULL, or contains any other value other than a valid UUID, then I want it to populate it with a UUID.  If there is already a valid UUID in the field, then I want it to do nothing.

I'm using ArcGIS Desktop 10.4.1 Standard License and I want to populate the field using the field calculator.  I'm also aware of GUID Fields in ArcGIS but they aren't suitable for the client's purposes.


texts: ogr2ogr, convert from shp to postgis error

texts:
I am trying to import a shapefile into PostGIS using the ogr2ogr tools that are built into the FWTools package.  In my import process I am receiving an error that I don't know how to deal with.

ERROR 4: Unable to open city.shx or city.SHX.
FAILURE: Unable to open datasource `city.shp' with the following drivers.

My code is:

ogr2ogr -overwrite -f "PostgreSQL" PG:"host=myhost user=myuser dbname=mydb password=mypass" city.shp

I have confirmed that the shapefile does not have the associated SHX file, but now what?  To make it worse, I have a whole lot of files to import that I assume will have the same issue because they also do not have an SHX.  Any thoughts on how I can resolve this problem?


texts: Determining if landlocked polygons in ArcGIS Desktop?

texts:
Imagine that I have the need to determine landlocked countries. Is there a way to do it automatically? Any known implementations?

The target platform on this will be ArcGIS.


texts:
I have been working with a set of data which mimics a shortest distance 'many to many' network analysis between COA population weighted centroids and education points in Bristol. My population centroids have been snapped to the road network, and I have tried various distance threshold within this section of code
-- ###
UPDATE bristol_coa11_centroids SET id = br.id
FROM bristol_roads_vertices_pgr br
WHERE ST_DWithin(bristol_coa11_centroids.geom, br.geom, 0.001);
-- ###
SELECT *
INTO test_route_v5
FROM pgr_dijkstra(
    'SELECT gid as id, source, target, length:: double precision AS cost
        FROM bristol_roads',
    ARRAY(SELECT id from bristol_coa11_centroids),
    ARRAY(SELECT id from bristol_education_v2), false) as be
    LEFT JOIN bristol_roads AS br ON be.edge = br.gid
    ORDER BY be.seq;

SELECT * FROM test_route_v5;

However, I have noticed that routes aren't being generated from all of the centroids (see image). Is there a way to get round this in Postgres or Arcmap first? I have measured distances from centroids and roads and determined different threshold against that but have seen no change. I need to expand the study size area so cannot afford to have a magnitude of points being missing. Thankyou!


texts:
I'm going to start to create programs to handle geographical data from a PostGIS server, and going to use GDAL.

According to your experience, with the best language to program with it? I can program in Perl and Python, but would like to know the language with more pieces of code to reuse, examples, libraries, etc.


texts: Checking for overlapping polygons in Oracle Spatial?

texts:
I have a table 1 which has a fields like store number and geometry and it just has one record in it. Now I want to check for this store number in table 2 and identify if there are any polygons overlapping to this particular store number polygon and insert them into the table 3.

Below is sql block I am trying to use.

Insert into table 3

select a.store_id,a.store_number,a.client_id,a.geometry
from  table 2 a,table 1 b
where a.client_id=b.client_id
and SDO_RELATE(a.Geometry, (select Geometry from table 1  where store_id
= 34746),'mask=anyinteract') = 'TRUE';

But it just returns 1 field in the new table for any store id I select.


texts: Geoserver publishes geography(Point) as Point, but geography(PointZ) as Geometry

texts:
CREATE TABLE schema.table
(
    location schema."geography(Point,4326)"(1107456) NOT NULL,
    location_z schema."geography(PointZ,4326)"(1107458)
)

I'm trying to figure out why Geoserver recognizes a geography(Point) as Point, but sees geography(PointZ) as Geometry.

I'm using Postgres 9.5, PostGIS 2.3, and Geoserver 2.11


texts: Geoserver Printing Module | MapFish print service won't print

texts:
I followed this documentation to use the printing service in geoserver.
I want to print the map layer in pdf format.

Below is my index.html for displaying the map and a button to print the map,

index.html:

 Export(pdf)

I have linked mapprint.js to apply the MapFish Print plugin.

mapprint.js

var specs={
 "layout":"A4 portrait",
 "srs":"EPSG:4326",
 "units":"degrees",
 "dpi":150,
 "mapTitle":"test",
 "layers":[
 {
    "baseURL":"http://localhost:9090/geoserver/BernHardt_Practical/wms",
    "opacity":1.0,
    "singleTile":false,
    "type":"WMS",
    "layers":["Nepal"],
    "format":"image/png"
 }
 ],
 "pages":[
 {
    "center":[84.2,28.2],
    "scale":2000000
 }
 ]
};

$("#export-pdf").click(function(){
var json = JSON.stringify(specs);
alert(json);
$.ajax(
{
  url:'http://localhost:9090/geoserver/pdf/print.pdf',
  type: 'POST',
  data:
  {
    spec: json
  },
  success: function()
  {
    alert("Success");
  },
  error: function()
  {
    alert("Failure");
  }
});
});

The map is displayed properly and ajax call is working fine but the problem is that it is not printing the pdf. I don't know where my mistake is.
According to the documentation, print.pdf command should have displayed the popup window for downloading the pdf file but nothing is happening.

My conf.yaml file is like this,

conf.yaml

# allowed DPIs
dpis: [75, 150, 300]

# the allowed scales
scales:
  - 25000
  - 50000
  - 100000
  - 200000

hosts:
  - !localMatch
    dummy: true

formats:
  - pdf
  - png

layouts:
  A4 portrait:
    mainPage:
      pageSize: A4
      rotation: true
      items:
        - !text
          text: '${mapTitle}'
          fontSize: 30
          spacingAfter: 30
        - !map
          spacingAfter: 30
          width: 400
          height: 400
        - !columns
          absoluteX: 410
          absoluteY: 310
          width: 100
          items:
            - !scalebar
              type: bar
              maxSize: 100
              barBgColor: white
              fontSize: 8
              align: right
        - !text
          font: Helvetica
          fontSize: 9
          align: right
          text: '1:${scale} ${now MM.dd.yyyy}'


texts: How to create hover effect on each polygon in mapbox map

texts:
function createPixelBox(jsonData){
    // a1=(x1,y2),a2=(x2,y2),a3=(x2,y1),a4=(x1,y1)

    var x,y,x1,x2,y1,y2;
    var geodata=[];
    var Length=jsonData.length;

          for(var i=0;i<Length;i++)
                {

                    x=jsonData[i].latitude;
                    y=jsonData[i].longitude;
                    console.log(x+"  "+y);

                    x1 = x - 0.0045; //lower x
                    x2=  x + 0.0045; //higher x
                    y1= y - 0.0045; //lower y
                    y2= y + 0.0045; //higher y

                   geodata.push({
                  "type" :"Feature",
                  "properties"  : {"description": jsonData[i].ndvi,
                      "pixel":jsonData[i]
                      },

                  "geometry"       : {
                     "type": "MultiPolygon",
                         "coordinates":[ [[[x1,y2], [x2,y2], [x2,y1], [x1,y1]]]]
                                    }
                               });

                   }

       prepareMap(geodata);

  }

      var valueSelector = 'ndvi';

         function   prepareMap(mapdata){

     L.mapbox.accessToken =  'pk.eyJ1IjoiZ2F1cmF2Y2F0c3RlY2giLCJhIjoiY2l1cGo0MTl4MDAxajJ1bng5a2xieTY0diJ9.0fDXNulI91U85ngSc4jGCg';
          var mapTwo = L.mapbox.map('map-two')
          .setView(defaultLocation,defaultZoom);
       L.mapbox.styleLayer('mapbox://styles/gauravcatstech/civ82vos7008f2is5x92cpvye').addTo(mapTwo);

  var myLayer = L.mapbox.featureLayer().addTo(mapTwo);

  var zoom = mapTwo.getZoom();
   console.log("Zoom level::::"+zoom);
  var states=mapdata;

    myLayer.setGeoJSON(states);

      myLayer.eachLayer(function(layer){

            console.log("Property in each layer");
            console.log(layer.feature.properties.pixel);

            var fillColor = "#B8B8B8"
            if (layer.feature.properties.pixel[''+valueSelector]=="38"){
              fillColor = "#ff0000";
            } else if (layer.feature.properties.pixel[''+valueSelector]<="250"){
              fillColor = "#ff3300";
            }else if ( layer.feature.properties.pixel[''+valueSelector]>="250"){
              fillColor = "#ffcc00";
            }

            layer.setStyle({
            fillColor: fillColor
                          });

                        });

          }

how would i apply hover effect on each polygon ....


texts: How to set the "LAYERS" parameter of the QGIS 3 algorithm "native:mergevectorlayers"

texts:
I do not understand how to set the LAYERS parameter for the QGIS 3 native: mergevectorlayers algorithm (for standalone python code).

Merge vector layers (native:mergevectorlayers)
...
----------------
Input parameters
----------------
LAYERS:  <QgsProcessingParameterMultipleLayers>
    Input layers

CRS:  <QgsProcessingParameterCrs>
    Destination CRS

OUTPUT:  <QgsProcessingParameterFeatureSink>
    Merged
----------------
Outputs
----------------
OUTPUT:  <QgsProcessingOutputVectorLayer>

The help processing.algorithmHelp('native:mergevectorlayers') says that this parameter has the data type QgsProcessingParameterMultipleLayers, but I do not understand how to use it.

Is this a list of the QgsVectorLayer objects?

Seems not because I can't execute the following code:

import sys
import qgis.core
qgis.core.QgsApplication.setPrefixPath(r'C:\OSGeo4W64\apps\qgis', True)
qgs = qgis.core.QgsApplication([], False)
qgs.initQgis()
sys.path.append(r'C:\OSGeo4W64\apps\qgis\python\plugins')
import processing
processing.core.Processing.Processing.initialize()
qgis.core.QgsApplication.processingRegistry().addProvider(qgis.analysis.QgsNativeAlgorithms())

input_files_list = [r"W:\MODIS\shp_files\MOD03_2018-01-01.shp",
                    r"W:\MODIS\shp_files\MOD03_2018-01-02.shp",]
layers_list = []
for num, file in enumerate(input_files_list):
    layer_name = 'layer_' + str(num)
    print(layer_name, file)
    layer = qgis.core.QgsVectorLayer(file, layer_name, 'ogr')
    print(layer)
    layers_list.append(layer)
    print('')

print('layers_list=', layers_list)
print('It seems that the script is stopped in the next line')
in_layers = qgis.core.QgsProcessingParameterMultipleLayers(layers_list)
print('This will not be printed')

output_file_path = r"W:\MODIS\shp_files\Merged.shp"
out_layer = qgis.core.QgsVectorLayer(output_file_path, 'out_layer', 'ogr')

params = {'LAYERS': in_layers,
          'OUTPUT': out_layer,}

feedback = qgis.core.QgsProcessingFeedback()

alg_name = 'native:mergevectorlayers'

print('='*50)
print(processing.algorithmHelp(alg_name))
print('='*50)
res = processing.run(alg_name, params, feedback=feedback)
print('='*50)
print(res['OUTPUT'])

qgs.exitQgis()

UPDATE 1:
I changed my stand-alone python code to check if I can pass parameters as a clear shapefile pathes, code execution gives this error:

Unable to execute algorithm
Incorrect parameter value for LAYERS

Here is new code:

import sys
import qgis.core
qgis.core.QgsApplication.setPrefixPath(r'C:\OSGeo4W64\apps\qgis', True)
qgs = qgis.core.QgsApplication([], False)
qgs.initQgis()
sys.path.append(r'C:\OSGeo4W64\apps\qgis\python\plugins')
import processing
processing.core.Processing.Processing.initialize()
qgis.core.QgsApplication.processingRegistry().addProvider(qgis.analysis.QgsNativeAlgorithms())

in_layers = [r"W:\MODIS\shp_files\MOD03_2018-01-01.shp",
             r"W:\MODIS\shp_files\MOD03_2018-01-02.shp",]
output_file_path = r"W:\MODIS\shp_files\Merged.shp"
params = {'LAYERS': in_layers,
          'OUTPUT': output_file_path,}
feedback = qgis.core.QgsProcessingFeedback()
alg_name = 'native:mergevectorlayers'
res = processing.run(alg_name, params, feedback=feedback)
print(res)
qgs.exitQgis()


texts: How long does it take to load the complete Tiger 2013 dataset into Postgis?

texts:
I'm loading the complete TIGER 2013 dataset into a Postgis database and finding that the process is taking an extraordinarily long time.  I was surprised by how long it took to load the TIGER 2012 dataset locally (over 60 hours), but the TIGER 2013 load is so long as to be almost unusable.

Loading a single state's worth of data (whether it be Alaska, Alabama, or Arkansas) is taking 12+ hours from my local machine on a high speed network into an Amazon RDS Postgres 9.3 instance with Postgis 2.1.  Note that this doesn't include downloading the files from ftp2.census.gov, since I pre-downloaded all the files.  This is all time taking to process the SHP files and insert them into Postgres.

I don't seem to be throttled by local resources - my CPU is nearly idle, and my load average is hovering around 2.  My network has plenty of spare bandwidth.

I can't imagine I'm throttled by network on the RDS instance, although it could be limited by CPU or memory (it's an m1.large RDS).  There is literally nothing else running against the database though, so the load should be minimal.

This situation is further complicated by the fact that the default TIGER scripts for Postgis don't lend themselves to parallelizability - they drop and create the tiger_staging schema, which makes it impossible to load states in parallel.  I'm planning to modify the scripts to allow parallel load, but before I do that I'd like to know if there's a more fundamental problem.

Has anyone done a full soup-to-nuts load of TIGER 2013 data?  If so, how long did it take?  Are there any configuration parameters I should be looking at that might be throttling this process?

Any advice or guidance would be appreciated.


texts: OpenGeo Suite 4.0.1 Activating SSL for remote connection

texts: I have installed OpenGeo Suite 4.0.1 on my Windows 7 OS and created a spatially database. Now i want to connect to my database remotely but then it seems like i have to generate a self-signed key. i read about SSL connection on http://suite.opengeo.org/opengeo-docs/dataadmin/pgDBAdmin/security.html but i have no idea where to write the codes for generating a self-signed key since my psql is not active anymore because i stoped the postgres service. If i write it in the command line it says openssl is not recognised as an internal or external command or batch file. Where can i write those codes to generate a self-signed key.


texts: WhiteboxGAT Map Gully Depth Parametrization

texts:
I'm trying to extract the ditch network from a DEM with the map gully depth tool in Whiteboxgat. Does anyone know about parameterization methodologies or where can I find them?
I found only two cases of application of this tool: Mahor et al. (2016) and Gustavsson and Selberg (2018)


texts: Is there a difference in results between Simplifying a geometry by different projections?

texts:
We have geometries in a PostGIS database projected in srid 4326.  We need to have these geometries simplified with ST_SimplifyPreserveTopology.

My question is whether this geometry needs to be transformed to a different projection--like 3400 in Alberta, Canada--to simplify it accurately by meters instead of degrees: if simplified by degrees, would it remove points more frequently on the parts of the geometry farther from the equator (because of the way degrees bend around the earth)?

Or would the function remove points in the same frequency and locations no matter the projection, if given the same region as a test?


texts: How to create circular polygon in openlayers 3?

texts:
So i would really like to modify this example:
http://openlayers.org/en/v3.0.0/examples/tissot.html?q=circle

Problem is that when I try to apply it to my map it doesn't work, probably because I use OSM style which is not spherical:

var map = new ol.Map({
  layers: [
    new ol.layer.Tile({
        source: new ol.source.OSM()
      }),
    new ol.layer.Vector({
      source: vectorSource
    })
  ],
  renderer: 'canvas',
  target: 'map',
  view: new ol.View({
    center: ol.proj.transform([2.1833, 41.3833], 'EPSG:4326', 'EPSG:3857'),
    zoom: 2
  })
});

And for the question: how to create circular polygon?
As I can see, there are two choices:

Somehow convert geom.Circle in geom.Polygon, which I can't do, based on my noob level
Create my own function to do that, something like openlayers 2:

OpenLayers.Geometry.Polygon.createRegularPolygon = function(origin, radius, sides, rotation) {
    var angle = Math.PI * ((1/sides) - (1/2));
    if(rotation) {
        angle += (rotation / 180) * Math.PI;
    }
    var rotatedAngle, x, y;
    var points = [];
    for(var i=0; i

Has anyone got to the same point in their GIS master life? Which one is better? or is there a hidden path, that only Wise and Strong ones can see? Help me please!


texts: Image Mosaic Layer Extents Changing Randomly

texts:
Any ideas what could cause the extents of a layer to change without me adjusting them?

Every so often I have to go back and reclick on the "compute from native bounds" under Lat/Lon Bounding Box:

The layer is an image mosaic layer in Geoserver. Nightly a script adds a new geotiff for the current day with the same extent as all previous tiffs inside of the image mosaic layer directory (it's a time series). I also add a new row to the postgis granule index table. But I copy the extend from the previous granules so they should all have the same bounding box. Would any of this possibly trigger Geoserver to change the Lat/Lon Bounding Box? Next time it happens I'll attach a screenshot where the numbers are off.

More Info - 5/10/2016

I have about 25 Geoserver layers. When I came in this morning, I checked the bounding box on all the layers and indeed the lat/lon bounding box on one of the layers had changed automatically. However it was a different layer than the four layers that this had happened to yesterday. The ones I fixed yesterday still look ok today. Here is a screenshot taken BEFORE I clicked "Compute from native bounds" to fix it.

and here is the screen shot after i clicked compute from native bounds

Also, I have this setup pretty much duplicated on a dev and ops server. Of course it only seems to be happening on the ops server :S


texts: ogr2ogr select feature class from multiple filegdb

texts:
From windows command shell how might I merge a single layer from multiple file-gdb's into another file-gdb using ogr2ogr? All input gdb's have identically named feature classes and attribute structures.

pseudo example:

ogr2ogr {input layername} D:\output.gdb x:\incoming\*.gdb

it's the layername part I'm having trouble with.


texts:
When exporting files from ArcPy anything that is created gets instantly imported in ArcMap so it can be used in following processes.

This is because ArcPy needs the files that uses imported to ArcGIS in order to find them.

How to import files to ArcGIS from a folder so I can use them in ArcPy without manually importing them?

For example a folder called files1

#set environment
arcpy.env.workspace = r"C:\Users\user\Desktop\axc\files1"
for i in arcpy.ListFeatureClasses():
    print(i)

This is for printing.
How can I import them to ArcGIS?


texts: writeOGR column limit error: "Creating Name field failed"

texts:
I have a SpatialPointsDataFrame with 1057 elements and 4760 columns and want to export it to a shapefile. When I try to do this with writeOGR() I get an error:

> writeOGR(pts,dsn="/home/joana/Desktop/",layer="PT", driver="ESRI Shapefile")
Error in writeOGR(pts, dsn = "/home/joana/Desktop/", layer = "PT", driver = "ESRI Shapefile") :
  Creating Name field failed

When I reduce the column number, everything works fine.But I really need all of the columns, is there any way to work around it?

I'm using Ubuntu 12.04/R/QGIS


texts: Looking for good load testing tools to test my web GIS

texts:
I am looking for good load testing tools to test my web GIS which runs on Apache/ESRI AU's Dekho application (Flash based, so would require AMF support within the load testing tool).

I am currently working SoapUI/LoadUI, but I am having a tough time creating steps for it, as there are numerous AMF calls from the client to the server. Some calls even have multiple calls embedded into it.

Has anyone used any load testing tool with GIS?


texts: How much time do I have untill the job status is no longer available?

texts: After launching a SubmitJob operation, how much time is left before the job status is no longer available?


texts: Clipping two rasters with several polygon features - saving output rasters according to a field name in the polygon layer?

texts:
I am an Economists and I am new to GIS. I am using ArcGIS for desktop 10.6. I am trying to clip a raster with several polygon features - saving output rasters according to a field name in the polygon layer.

I used the code by Midavalo. The code is running perfectly but the problem is that it gives me an output of light raster for one province out of 200 I have in my field, it looks like it does not loop through i. It has been 2 months working on this but I get no help.

import arcpy, os

arcpy.env.overwriteOutput = True

arcpy.env.workspace = "C:\Users\Somaqaqa\Desktop\jupyter_python-tutorial1\installPycharm_pythonscript+\creating a file database\GiSPython_China\just_rastes"

clipshapefile = r"C:\Users\Somaqaqa\Desktop\jupyter_python-tutorial1\installPycharm_pythonscript+\creating a file database\GiSPython_China\Data\chn_admbnda_adm1_ocha.shp"

rasterlist = arcpy.ListRasters()

clipshapes = set() with arcpy.da.SearchCursor(clipshapefile, ['ADM1_EN']) as cursor:

for row in cursor:
     clipshapes.add(row[0])
Loop through input rasters, and clip by each shape in the input shapefile
for raster in rasterlist:

rastername, rasterext = os.path.splitext(raster)
for i in clipshapes:

    newRaster = "{}_FID{}.tif".format(rastername, i)
    newRasterPath = os.path.join(outputWorkspace, newRaster)
    Cliplayer = "Ras {}".format(raster)

 if arcpy.Exists('ClipLayer'):

arcpy.Delete_management('ClipLayer')
arcpy.MakeFeatureLayer_management(clipshapefile, Cliplayer, " ADM1_EN = '{}'".format(i))

arcpy.Clip_management(raster, "#", newRasterPath, Cliplayer, "0", "ClippingGeometry",
"MAINTAIN_EXTENT")
print raster

print i


texts: Set Raster Cells Based on NoData from other Raster using ArcGIS Desktop?

texts:
I have just no clue on this problem. I want one raster to set its values to 0 or NoData based upon another raster noData or 0 values. It's like a reclassify, but between two different raster's. I have a picture here to help illustrate:

I have two raster here: one is a depth raster (the one with the reds and has a flow look to it) and the other is a velocity raster (the one with the straight breaks).

I'd like to have the Velocity raster only exist where the Depth raster exists. I think this needs to be done in raster calculator or a clip, but so far I haven't found any way to do this.


texts:
I made a surface plot with matplotlib, mplot3d and gdal. Here is the code:

import gdal
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np

# maido is the name of a mountain
# tipe is the name of a french school project

# 1) opening maido geotiff as an array
maido = gdal.Open('dem_maido_tipe.tif')
dem_maido = maido.ReadAsArray()

# 2) transformation of coordinates
columns = maido.RasterXSize
rows = maido.RasterYSize
gt = maido.GetGeoTransform()

x = (columns * gt[1]) + gt[0]
y = (rows * gt[5]) + gt[3]

X = np.arange(gt[0], x, gt[1])
Y = np.arange(gt[3], y, gt[5])

# 3) creation of a simple grid without interpolation
X, Y = np.meshgrid(X, Y)

# 4) deleting the "no data" values

# delete the last column
dem_maido = np.delete(dem_maido, len(dem_maido)-1, axis = 0)
X = np.delete(X, len(X)-1, axis = 0)
Y = np.delete(Y, len(Y)-1, axis = 0)

# delete the last row
dem_maido = np.delete(dem_maido, len(dem_maido[0])-1, axis = 1)
X = np.delete(X, len(X[0])-1, axis = 1)
Y = np.delete(Y, len(Y[0])-1, axis = 1)

# 5) plot the raster
fig, axes = plt.subplots(subplot_kw={'projection': '3d'})
surf = axes.plot_surface(X, Y, dem_maido, rstride=1, cstride=1, cmap=cm.gist_earth,linewidth=0, antialiased=False)

plt.colorbar(surf)  # adding the colobar on the right

plt.show()

Here is the GeoTiff file that I used (it is uploaded on my own google account don't worry):

https://drive.google.com/file/d/0B7P95aWmH4DUQk9SbzhNUVNINGs/view

And here is what I get:

I am happy with the result but when I try to move it, it takes to much time.

What can I do to make it faster ?


texts: Is there any web-gl js libraries for creating 3D bar overlap city map?

texts:
Is there any some web-gl js libraries for creating 3D bar overlap city map?

More details:
http://ny.curbed.com/archives/2015/01/09/3d_map_shows_where_nyc_housing_is_priciest_cheapest.php

https://medium.com/re-form/nycs-housing-cost-myth-9dce6052c139

EDIT:
To narrow the answer set, an interactive web mapping solution without data uploading will be the best answer.


texts: Filter ImageCollection by year and create one mean Image for each year

texts:
I'm trying to calculate the mean NDVI for a polygon in an urban area.
I have a code that takes Landsat7 images, clip them according to the geometry and knows to select only images with certain number of pixels.

The next step is that I want to take those images that are ranging from 1999-2013, and for each year to create one mean image. so in the end instead of for example 67 images, I'll have one image per year.
 I have faced two problems:

I have tried to use :

var img2000=withNDVI.filter(ee.Filter.calendarRange,(2000,2001,'year')).mean();

in order to get one mean image from 2000, but I get the next error message:

  Received too many arguments to function filter(). Expected at most 1
  but got 2.

and I'm also not sure that if I choose 2000-2001, if it knows that I want only 2000 or it will take also 2001.

I'm lazy. I don't want to write this line so many times and I want to have function that does it -that knows to calculate mean for each year but not sure how to do it, should I start with defining a variable for the function?

Here is the code I have:

var geometry=Osaka;
Map.centerObject(geometry,12);
//create Landsat7 database
/**
 * Function to mask clouds based on the pixel_qa band of Landsat SR data.
 * @param {ee.Image} image Input Landsat SR image
 * @return {ee.Image} Cloudmasked Landsat image
 */
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2).divide(10000)
  .copyProperties(image, ['system:time_start']);
};

//Create LANDSAT7 dataset

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filterDate('1999-01-01', '2013-04-29')
                  .select('B3','B4','pixel_qa')
                  .filterBounds(geometry)
                  .map(cloudMaskL457);

//clip the dataset according to the geometry
var clippedCol=dataset.map(function(im){
  return im.clip(geometry);
});

// Get the number of images.
var count = clippedCol.size();
print('Count: ',count);

//function to calculate NDVI in LANDSAT7

var addNDVI = function(image) {
  var NDVI = image.normalizedDifference(['B4', 'B3'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(NDVI);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var NDVIcolor = {
  min: 0,
  max:1,
  palette: ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'],
};

var img2000=withNDVI.filter(ee.Filter.calendarRange,(2000,2001,'year')).mean();


texts:
I am using ArcGIS Pro 2.6.1.
I have massive map of roads in multiple road files. I need to extend lines to their intersect point. Intersect tool does not seem to work to my satisfaction. Any other suggestions.
Should I combine multiple files into one?
I have a layer of points that can be extended to, but could not find anything.


texts: Producing raster A + B where A and B intersect, and A where they do not, using QGIS?

texts:
I have a raster A (a heightmap created from a contour file). I also have created a vector layer with some polygons which I rasterized into raster B. Raster B is completely contained within the extents of raster A.

Raster A has values at each pixel. Raster B, on the other hand, has a bunch of empty area (which I believe is represented as 0 or nodata).

I'd like to produce raster C, such that C = A + B, where A and B intersect, and C = A where A and B don't intersect.

The problem I'm running into is that when I add the two rasters together using the raster calculator, the result only has values where the two rasters intersect, but show nothing where they do not intersect. I'm guessing this is due to how QGIS treats nodata values.

How can I achieve what I'm trying to do?


texts: Number points along polyline while keeping information on preceding point

texts:
I'm working with river networks and need to create ID for points along it. Trick is that points need to have its own ID and also keep the ID of the preceding point (in  two separate columns of course).

It would look like this:

Are there any tools/plugins in QGIS that I could use to create such data? Preferably without coding.


texts: OpenLayers canvas renderer - redraws all features on selecting a single feature

texts:
I am currently trying to learn more about OpenLayers, and I'm coming up against performance / behaviour issues that I wouldn't expect to see. To elaborate...

I have combined this example for generating lots of random points with this example on highlighting / changing the icons of features on selection.

The result is available in this JSFiddle.

The (very basic) application prioritises the Canvas renderer, so that the points layer is defined like this:

var layer = new OpenLayers.Layer.Vector('points', {renderers: ['Canvas', 'SVG', 'VML'],
styleMap: new OpenLayers.StyleMap({
    "default": new OpenLayers.Style(OpenLayers.Util.applyDefaults({
        externalGraphic: "img/marker-green.png",
        graphicOpacity: 1,
        rotation: 0,
        pointRadius: 10
    }, OpenLayers.Feature.Vector.style["default"])),
    "select": new OpenLayers.Style({
        externalGraphic: "img/marker-blue.png"
    })
})});

The Canvas renderer is nice because zooming and panning works well with thousands of points (unlike SVG) but clicking / feature selection is much slower than SVG. I was surprised to see that all feature icons are redrawn before the selected feature's icon changes (to marker-blue.png). This is especially apparent in Firefox, but less noticeable in Chrome and too fast to see in IE9.

This question on SO suggests that there could be an issue with how OpenLayers is interacting with the Canvas, but unfortunately I don't know enough about HTML5 or OpenLayers (yet) to be sure.

Does anyone have experience with OpenLayers / Canvas renderer? I am interested in any ways of speeding up feature selection when there are lots of features - whether it's the time taken to figure out which feature was selected, or simply changing the selected feature's icon.

Any help much appreciated.


texts: COLLADA files and geographic coordinates

texts:
I am failing to understand something here, from the arcgis help page:

  "The Multipatch To COLLADA tool can be used to export multipatches to
  COLLADA files on disk. COLLADA is a 3D model format that is supported
  in many 3D applications, including Google Earth. If the multipatch
  feature class you are converting is in a projected coordinate system,
  the tool will create a KML file in the output folder which will assign geographic coordinates to the 3D models in the COLLADA files"

When I convert a multipatch to COLLADA i see nothing more than the .DAE files. Anyone encounter this before?


texts: Overlaying thousands of points of a route from shapefile and dbf on google maps

texts:
I have to overlay thousands of points(data taken every second) of a route from shapefile and related data for each points in dbf. I have to show point related data in a popup on click of a point. I will be getting few shapefiles every day and overlay data of few of the shapefile on the map as per request. Over 2-3 years I will have multiple thousand shapefiles. I have tried two different approaches as follows.

Parse the shapefile with PHP classes at request time with pointer of shapefile in database table (postgres db) and pass all the data to browser in JSON. It becomes very slow in overlaying points on google maps as I zoom in on to the route.

Another approach I used is storing each points(long/lat) with corresponding data from dbf as separate rows and connect my db(postgres) to geoserver. Using geoserver to serve data tile-by-tile. Making request to geoserver to show cached point data on click of a point, but then I am worried about dealing with millions of rows (thousands of points of a shapefile with multiple thousand shapefile over couple of years) in postgres table.

Do you have any suggestions for some other approach or improving one of the current approaches?


texts: Converting OpenStreetMap data to shapefile using OGR?

texts:
How would I reliably convert OSM-data (.osm, .kml, .gpx or .geojson) to .shp?
ogr2ogr seems have issues:

geojson:

$ogr2ogr -f "ESRI Shapefile" destination_data.shp export.geojson
ERROR 6: Can't create field alt_name31 in Shape DBF file. Maximum record length reached.
Warning 6: Normalized/laundered field name: 'alt_name:fa' to 'alt_name31'

kml, gpx:

$ogr2ogr -f "ESRI Shapefile" destination_data.shp export.kml  -lco ENCODING=UTF-8
Warning 6: Normalized/laundered field name: 'Description' to 'Descriptio'
Warning 1: Value '...' of field Descriptio has been truncated to 254 characters.

And encoding is off.


texts: Differences in raster extents in gdalinfo, SAGA & Quantum after resampling

texts:
I recently resampled a raster file using SAGA's Resampling module. The aim was to widen the extent of the new raster to two decimal places, where the original raster has three decimal places. The resolution of the raster stays the same with 0.02. I take the extent of the original raster and calculate the extent of the additionally needed input-raster from it.

For example - Orginal raster upper left corner:

This is what gdalinfo tells me:

589262.561, 5955004.195

What  I see with SAGA hard zoom to displayed raster:

589262.561, 5955004.195

What  I see in SAGA-Metadata instead:

589262.571, 5955004.185

What  I see with Quantum hard zoom to displayed raster:

589262.561, 5955004.195

What  I see in Quantum Metadata:

589262.5609, 5955004.1949

The given numbers of lines and columns are all the same from gdalinfo to SAGA and Quantum. Why is there such a difference in displaying the raster in SAGA?

After resampling with the calculated input grid I get the following:

gdalinfo:         589262.550, 5955004.210
SAGA-zoom:        589262.570., 5955004.190
SAGA-metadata:    589262.56, 5955004.2
Quantum-zoom:     589262.570., 5955004.190
Quantum-metadata: 589262.550, ,5955004.209

Again such differences between displaying the raster and its actual metadata. However, the amounts of lines and columns are identical again, though the other corners also differ slightly.

Is there generally such inaccuracy in these programs? Is it somehow caused by the resampling module? Which information should I believe for further use of the raster?

A lot of information and questions--I would be grateful to get an answer anyway.


texts: ArcPy "Cannot open table for Load" error

texts:
I have been using some slightly modified code based on socalgis.org to extra data from map services. I have had no problems until this week where I am now getting the error

  RuntimeError - Cannot open table for Load

I am getting this on services that had worked fine in the past. For example the following URL works fine with no error whilst this fails with "Cannot open table for Load" error. Both are "esriGeometryPolygon" with a small number of feature on the same server.

Record extract limit: 1000
Number of target records: 403
Gathering records...
  OBJECTID >= 1 and OBJECTID <= 403
  Query: https://services.gis.ca.gov/arcgis/rest/services/Environment/Weather_stations/MapServer/2/query?where=OBJECTID >= 1 and OBJECTID <= 403&returnGeometry=true&outFields=*&f=json

Traceback (most recent call last):

  File "<ipython-input-39-454949d8d113>", line 1, in <module>
    runfile('[Path to code]', wdir='[Path to code]')

  File "C:\Users\User\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone11\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 827, in runfile
    execfile(filename, namespace)

  File "C:\Users\User\AppData\Local\ESRI\conda\envs\arcgispro-py3-clone11\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 110, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "[Path to code]", line 50, in <module>
    fs[i].load(urlstring)

  File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\arcobjects\arcobjects.py", line 421, in load
    return convertArcObjectToPythonObject(self._arc_object.Load(*gp_fixargs(args)))

RuntimeError: RecordSetObject: Cannot open table for Load

The error seems to resolve around fs[i].load(urlstring). My code below:

import arcpy
import urllib
import json

arcpy.env.overwriteOutput = True
baseURL = "[URL to map service]"
fields = "*"
outdata = "[path\to\gdb]"

# Get record extract limit
urlstring = baseURL + "?f=json"
j = urllib.request.urlopen(urlstring)
js = json.load(j)
maxrc = int(js["maxRecordCount"])
print (("Record extract limit: %s" % maxrc))

# Get object ids of features
where = "1=1"
urlstring = baseURL + "/query?where={}&returnIdsOnly=true&f=json".format(where)
j = urllib.request.urlopen(urlstring)
js = json.load(j)
idfield = js["objectIdFieldName"]
idlist = js["objectIds"]
idlist.sort()
numrec = len(idlist)
print (("Number of target records: %s" % numrec))

# Gather features
print (("Gathering records..."))
fs = dict()
for i in range(0, numrec, maxrc):
  torec = i + (maxrc - 1)
  if torec > numrec:
    torec = numrec - 1
  fromid = idlist[i]
  toid = idlist[torec]
  where = "{} >= {} and {} <= {}".format(idfield, fromid, idfield, toid)
  print (("  {}".format(where)))
  urlstring = baseURL + "/query?where={}&returnGeometry=true&outFields={}&f=json".format(where,fields)
  print (("  Query: %s" %urlstring))
  fs[i] = arcpy.FeatureSet()
  fs[i].load(urlstring)

# Save features
print (("Saving features..."))
fslist = []
for key,value in fs.items():
  fslist.append(value)
arcpy.Merge_management(fslist, outdata)
print (("Done!"))

I am using Python 3.6 and Spyder installed in a ArcGIS Pro environment. How can I  resolve this?

UPDATE
The service with the error seems to be caused by an incorrect query - the above example being generated by URLSTRING

https://services.gis.ca.gov/arcgis/rest/services/Environment/Weather_stations/MapServer/2/query?where=OBJECTID >= 1 and OBJECTID <= 403&returnGeometry=true&outFields=*&f=json

Which gives

{"error":{"code":400,"message":"Failed to execute query.","details":[]}}


texts: Resample binary raster to give *proportion* within new cell window

texts:
I would like to go from a binary raster of forest/non-forest at 30m resolution to a raster at 240m resolution with the value of each cell being the proportion of that cell forested - i.e. the data values for the 240m cells would go from 0 (where all cells in the 30m raster were zero/non-forest) through 0.5 (where half the 30m cells were forest, half were non-forest) to 1 (where all cells in the 30m raster were forested).

Although it is has been suggested not to use bilinear interpolation on discrete data, as far as I can tell the result with binary 0/1 data would be to give an average value (i.e. a proportion between 1 and 0). Is this a sensible way to do it, or is there a better method?

I can use Arc, QGIS and Idrisi.


texts: Deleting file geodatabase using ArcPy with ArcGIS Pro?

texts:
I am trying to delete a geodatabase using arcpy but I don't know how to do it. I have tried with 'Delete_management' but it doesn't work. Any ideas?

path = D:\ArcGIS_Pro_Projects\KML files\geodatabase_name.gdb
arcpy.Delete_management(path)


texts: Count touching and near points for every polygon (QGIS)

texts:
I have a polygon and a point layer. The polygons contains buildings, the points should be the edges of the buildings. A building is "good" as much as points are touching or in the near (buffer 0.2 meter) of the building. There are no other references between the layers.

I need the number of touching or intersecting points and points within a near of 20cm. With that result I'm able to say, which building have got the best quality...

The result could be written into a virtual field of the building-polygon.

I know, SQL could do the thing easily. But I need a QGIS (aggregate?) way because the input data will be provided as shapefile or GeoJSON.


texts: QGIS field header - change default number of digits

texts: In QGIS the attribute table of an imported CSV has some field headers with more than 9 digits. In a derived shapefile, field headers, where source header is more than 9 digits, get truncated. Because the headers are different when the source is updated the fields in the shapefile with truncated headers are not populated (they are null). Is there a recommended way to deal with this?


texts: Adding data to subgroup layers with same name in ArcPy?

texts:
Assuming that I have an idealized .mxd TOC structure with group layers that contain group and subgroup layers which have the same name like in the following example:

Layers

GroupLayer1

SubGroupLayer1

Raster1.tif
Raster33.tif

SubGroupLayer2

Raster44.tif

GroupLayer2

SubGroupLayer1

Raster2.tif
Raster4.tif

SubGroupLayer2

Raster7.tif

Is it possible to add data using an 'absolute' layer path with ArcPy?

For example, if I want to add a raster layer SampleRaster.tif to the following place GroupLayer2/SubGroupLayer2 : how could I do this in ArcPy?

I know how to add a group layer and then add the specific layer (e.g. solution in this question Adding new group layer with ArcPy?) but have not found a solution how to add it if two layers have the same name.

My code that I have so far is pasted below

 # add rasters to .mxd project
    import arcpy

    arcpy.env.workspace = r'D:\Data\Rasters'
    rasters = arcpy.ListRasters("*", "TIF")

    for raster in rasters:
        # add raster to .mxd
        mxd = arcpy.mapping.MapDocument(r"D:\Data\Test.mxd")
        df = arcpy.mapping.ListDataFrames(mxd)[0]

        # here I need to list the layers in the subgroup "SubGroupLayers2", but what to do when the subgroup layers names are equal? Is there a way to indicate an absolute path (e.g. Layers/GroupLayer2/SubGroupLayer2 here?)
        targetGroupLayer = arcpy.mapping.ListLayers(mxd, "SubGroupLayer2", df)[0]

        # I am then using another layer to apply the symbology to the raster and update the layer
        newlayer = arcpy.mapping.Layer(raster)
        newlayer.visible = False
        arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, newlayer, "BOTTOM")
        updateLayer = arcpy.mapping.ListLayers(mxd, raster, df)[0]
        sourceLayer = arcpy.mapping.Layer(r"D:\Data\Colormap.lyr")
        arcpy.mapping.UpdateLayer(df, updateLayer, sourceLayer, True)
        mxd.save()


texts: How to use buffers to estimate the attribute data for a given polygon?

texts:
I have two layers: point and a polygon.

The point layer has lat/long and I have created buffers of a 1/4 mile.

The polygon layer has data such as age and income.

Is there a plug-in that can interpolate the attribute data (e.g. age or income) contained within a buffer area?

I recall this being a plug-in in ArcView, but I have no clue of its contemporary (or the process, limitations) in QGIS.


texts:
I have added labels showing the District names of Lesotho, but the highlighted feature is not showing a label until I zoom in or out.

Zoom to country extent and label is missing:

The label shows when zoomed in:

"Show all labels for this layer" has been checked.

Is this a bug with version 3? I am currently using 3.0.3


texts:
I've been looking into tools for displaying dynamic terrain (x,y,z(t)) via the browser and I just recently came across cesium.js which looks like a wonderful tool at first glance. Many of the demos show static terrain which is great, but the data sets I'll be working with are all dynamic. My original plan was to do this with three.js and d3.js but after looking at cesium.js, I would say, "I'd be reinventing the wheel" (which would likely come out square).

I read up on the quantized-mesh and it sounds like the mesh is optimized based on the elevation. While that results in an optimized user experience for static displays, changing the height would require a recalculation of the mesh, which is practically of no use for continuous dynamic visualizations.

If  cesium.js does not support dynamic terrain of the form (x,y,z(t)) out of the box, what processes would I have to perform to display it through cesium.js?


texts: Generating coordinates with attributes from a georeferenced map

texts:
I've got a (physical) hydrographic map that I'm intending to scan and georeference into QGIS. On this map are hundreds of point depths measurements and what I want is to obtain an  array with the coordinates and depth for each point.

Is there a straightforward way of doing this in QGIS or ArcGIS (e.g. clicking on a point to copy the coordinates into an attribute table where I can manually enter the depth for that point)?


texts: How would I configure a geoserver store using jTDS?

texts:
I'm trying to create a GeoServer docker container, which accesses a SQL Server database, and would prefer to use the jTDS driver over the Microsoft one--besides the fact that it's supposed to be faster, installing the Microsoft driver into a docker container is a pain (the URL keeps changing).

My GeoServer install includes the sqlserver plugin, and so consequently, among Stores I see options for Microsoft SQL Server (JNDI) and Microsoft SQL Server (JTDS Driver) (JNDI).

When I add libtds-java to my (ubuntu-based) docker, and copy the jtds.jar file to my $GEOSERVER_HOME/webapps/geoserver/WEB-INF/lib/ directory, then I get a Microsoft SQL Server (JTDS Driver) option, so that would appear to have found the appropriate driver.

But any attempt to create a Store results in the error Error creating data store, check the parameters. Error message: Could not acquire data access. I see this error mentioned for many reasons on the mailing list, and it appears usually to be a permissions problem on the data directory, but the docker is running as user www-data, everything in the data directory is owned by www-data, and if I exec bash in the docker I can manually create files in the data directory. Or I can run the docker as root, and still see the same error.

What's more, if I install the Microsoft driver, I can create stores using it.

Yet, I've seen mention that the jTDS driver works, so what do I need to do to configure it?

Tomcat solutions are acceptable. I'm not totally tied to the Jetty implementation. I suspect, in the end, adding the GeoServer war to a standard Tomcat docker would be simpler, anyway (though I'm not sure about the plugins).

And yeah, the real solution would be to use PostGIS. Why we use SQLServer is beyond me.


texts: Strange output from exif tiff data, degrees decimal mins

texts:
I am extracting GPS data from .tif images using pythons exifread. It reports the following:

After a bit of jiggery-pokery I get the following:

6 47.04774 W 49 13.08169 N

Am I correct in assuming that the division of 1308169/100000 is a correct method - it does match the windows properties which list the GPS as deg,min,sec?

And, how might I convert this format to decimal degrees?

Reproducible in Python:

print(long_num)
print(long_ref)
print(lat_num)
print(lat_ref)

gives:

runfile('blah/blah/blah', wdir='blah/balh')
[6, 2352387/50000, 0]
W
[49, 1308169/100000, 0]
N

obtained from:

import exifread

f = open(test_file, 'rb')

tags = exifread.process_file(f)

long_num  = str(tags.get('GPS GPSLongitude'))
long_ref  = str(tags.get('GPS GPSLongitudeRef'))

lat_num  = str(tags.get('GPS GPSLatitude'))
lat_ref  = str(tags.get('GPS GPSLatitudeRef'))


texts: Extracting land cover type in Google Earth Engine?

texts:
I have a question related to Drawing polygon and extracting land cover inside using Google Earth Engine?

I'm trying to work out how to calculate habitat type area from within a specific polygon in google earth engine using Copernicus data. I can do both of these tasks, just not together.

var s1 = ee.ImageCollection("COPERNICUS/S1_GRD"),
S2 = ee.ImageCollection("COPERNICUS/S2"),
fc = ee.FeatureCollection("ft:1w_ciwxR9FKm0rzMH8oJ3MdyE-1U15PfiOFlGW5Ib"),
fc2 = ee.FeatureCollection("ft:18Xvf7jG2EI4mJMwj0oLV2u-05bAPiwJDaxSRy2ep");
//Set a point of interest
var aoi = fc;
//Set bands of interest
var bands = ['B2', 'B3', 'B4','B8', 'B11', 'B12'];

//Import and filter Sentinel 2 collection
var S2 = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2016-07-01', '2016-08-01')
.filterBounds(aoi)
print(S2);
var image = S2.mosaic();

Map.addLayer(S2, {'bands': 'B12,B8,B4','min': 0,'max': 4000}, 'Sentinel 2A         SWIR,NIT,R', 0);
var training = image.sampleRegions(fc2, ["info"], 30);
var classifier = ee.Classifier.cart()
.train(training, "info", image.bandNames());
var out = image.classify(classifier);
Map.addLayer(out, {palette: [
'ff0000',// cleared red
'ffcc66',// cultivated
'387242',// forest
'99ff33',// grass
'e0e0d1',// rock
'e0e0d1'// rock grass
], min:1, max:6});

// Create a geometry representing an export region.
var geometry =         ee.Geometry.Rectangle([36.977296602,-15.364411987,37.057018981,-15.447436229]);

Map.addLayer(ee.Image().paint(fc, 1, 3));
var fc = fc2;
Map.addLayer(fc2);

//Create a band for each veg type
var class1 = out.eq(1).select([0],['cleared']);
var class2 = out.eq(2).select([0],['cultivated']);
var class3 = out.eq(3).select([0],['forest']);
var class4 = out.eq(4).select([0],['grass']);
var class5 = out.eq(5).select([0],['rock']);
var class6 = out.eq(6).select([0],['rockgrass']);
// Concatanate the bands into a single image
var total = ee.Image.cat([class1,class2,class3,class4,class5,class6]);

// Extract the landcover band
var landcover = image.select('landcover');

 // Clip the image to the polygon geometry
 var landcover_roi = landcover.clip(geometry);

 //Calculate the total area of each veg type in meters^2
 var count = total.multiply(ee.Image.pixelArea());
 print(total);
 print(count);

 var area = count.reduceRegion(ee.Reducer.sum(), geometry, 30, null, null,        false, 1e10, 1)
print(area);
var cleared = ee.Number(area.get('cleared'));

Export.image.toDrive({
image: out,
description: 'imageToDriveDeforestation',
scale: 30,
region: geometry
 });