texts:
Why is my map not in the center when I load my map in OpenLayers? And when I drag my map, it doesn't move sometimes.

var maxExtent = new OpenLayers.Bounds(103.30029292030527, 24.39151279999997, 109.88812867358875, 29.45420120000002);
        function init() {
            var jsonp = new OpenLayers.Protocol.Script();
            jsonp.createRequest(layerURL, {
                f: 'json',
                pretty: 'true'
            }, initMap);
        }

        function initMap(layerInfo){
            var baseLayer = new OpenLayers.Layer.ArcGISCache("AGSCache", layerURL, {
                isBaseLayer:true,
                layerInfo: layerInfo,
                StartBounds: maxExtent,
                transitionEffect:'resize'

            });
            /*
             * Make sure our baselayer and our map are synced up
             */
            console.log( baseLayer.maxExtent);
            map = new OpenLayers.Map('map', {
               // restrictedExtent: maxExtent,
                //startBounds: maxExtent,
               startBounds: maxExtent,
                maxExtent: baseLayer.maxExtent,
                units:"degrees",
                resolutions: baseLayer.resolutions,
                numZoomLevels: baseLayer.numZoomLevels,
                tileSize: baseLayer.tileSize,
                displayProjection: baseLayer.displayProjection
            });
            map.addLayers([baseLayer]);

            /*
            *  delete the default controls
            * */
            var zoomControl = map.getControlsByClass('OpenLayers.Control.Zoom')[0];
            map.removeControl(zoomControl);
            map.addControl(new OpenLayers.Control.MousePosition());
            map.addControl(new OpenLayers.Control.Attribution());
            map.addControl(new OpenLayers.Control.Navigation({dragPanOptions: {enableKinetic: true}}));
            map.addControl(new OpenLayers.Control.PanZoomBar({zoomWorldIcon:true}));
          // map.addControl(new OpenLayers.Control.LoadingPanel());
            if (!map.getCenter()) {
                map.zoomToMaxExtent();
            }
            /*zoom the map to a special level
             map.zoomTo(3,(107.4833233,26.23578279));
            map.setCenter(new OpenLayers.LonLat(107.4833233, 26.23578279), 2);*/
             map.zoomTo(0);
            map.addControl(new OpenLayers.Control.LayerSwitcher());
            map.addControl(new OpenLayers.Control.MousePosition() );
            //map.zoomToExtent(new OpenLayers.Bounds(-8341644, 4711236, -8339198, 4712459));
          //  map.zoomToExtent(new OpenLayers.Bounds(103.30029292030527, 24.39151279999997, 109.88812867358875, 29.45420120000002));
        }


texts: How to draw one rectangle in ArcObjects?

texts:
Is there button in ArcObjects "draw rectangle" such as in openlayers?

Update

Well, I use this code for drawing:

ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = Map1.ActiveView.ScreenDisplay;
screenDisplay.StartDrawing(screenDisplay.hDC, (System.Int16)
    ESRI.ArcGIS.Display.esriScreenCache.esriNoScreenCache);
ESRI.ArcGIS.Display.IRgbColor rgbColor = new ESRI.ArcGIS.Display.RgbColor();
rgbColor.Red = 255;

ESRI.ArcGIS.Display.IColor color = rgbColor;
ESRI.ArcGIS.Display.ISimpleFillSymbol simpleFillSymbol = new
    ESRI.ArcGIS.Display.SimpleFillSymbol();
simpleFillSymbol.Color = color;

ESRI.ArcGIS.Display.ISymbol symbol = simpleFillSymbol as
    ESRI.ArcGIS.Display.ISymbol;
ESRI.ArcGIS.Display.IRubberBand rubberBand = new
    ESRI.ArcGIS.Display.RubberEnvelope();
ESRI.ArcGIS.Geometry.IGeometry geometry = rubberBand.TrackNew(screenDisplay,
    symbol);
screenDisplay.SetSymbol(symbol);
screenDisplay.DrawRectangle(geometry as ESRI.ArcGIS.Geometry.IEnvelope);
screenDisplay.FinishDrawing();
...
ArcGisToolbar.CurrentTool.Deactivate();
ArcGisToolbar.CurrentTool = null;
ICommand command = this.ArcGisToolbar.CommandPool.get_Command(3); // button "pan"
Map1.CurrentTool = (ITool)command;

But there is a problem. When I activate this instrument (button "select") again, I already start selecting the area (from the point, where I begin drawing the polygon)


texts: Updating the file path for a Google Earth layer

texts:
I have exported a lot of our data from ArcMap to Google Earth. Most of our company prefers to use Google Earth for basic searches. The users within the company have downloaded a "Master" Google Earth file which has all of our GIS layers.

I have updated some of the layers in ArcMap and have exported them as KMZ. If I change the file path of the layers to show updated data in the "Master" on my PC, will all users in the company automatically get the updates the next time they start Google Earth?


texts: Difference between Priority in the Positioning dialog and Z-Index in the Rendering dialog in QGIS

texts:
I studying the layer labeling options in QGIS and encountered two options which seem to be the same. If I understand the QGIS documentation correctly, the option "priority" within "Positioning" determines which label to render when there are two Labels overlaping.

The Z-Index Option in "Rendering" seems to be doing the same.

What is the difference between these two options?


texts:
I am trying to set up a tile map service on Apache. I tried following this:
Setting Up a TMS Tile Map Service on a Web Server but when I go to load the page it's not loading the images in my /var/www/TMS/1.0.0/ne1 directory returning "Failed to load resource: the server responded with a status of 404 (Not Found) http://mysite-name.com/TMS/1.0.0/ne1/2/0/1.png" in the error console.

In the Apache config file:

Alias /map/ "/var/www/TMS/1.0.0/ne1/openlayers.html"
<Directory /var/www/TMS>
    Options Indexes MultiViews FollowSymLinks
    AllowOverride None
    Order Allow,Deny
    Allow from all
</Directory>

root.xml file:

<?xml version="1.0" ?>

<Services>

  <TileMapService title=" Tile Map Service" version="1.0.0" href="http://mysite-name.com/TMS/1.0.0/service.xml" />

</Services>

service.xml file:

<?xml version="1.0"?>

<TileMapService version="1.0.0" services="http://mysite-name.com/TMS/root.xml">

<TileMaps>

   <TileMap srs="EPSG:4326" title="ne1" profile="geodetic" href="http://mysite-name.com/TMS/1.0.0/ne1/tilemapresource.xml"/>

</TileMaps>

</TileMapService>

tilemapresource.xml:

<?xml version="1.0" encoding="utf-8"?>
    <TileMap version="1.0.0" tilemapservice="http://mysite-name.com/TMS/1.0.0/service.xml">
      <Title>test.vrt</Title>
      <Abstract></Abstract>
      <SRS>EPSG:4326</SRS>
      <BoundingBox minx="-125.00000000000000" miny="20.00330020777438" maxx="-54.99858571732217" maxy="50.00000000000000"/>
      <Origin x="-125.00000000000000" y="20.00330020777438"/>
      <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
      <TileSets profile="geodetic">
        <TileSet href="2" units-per-pixel="0.17578125000000" order="2"/>
        <TileSet href="3" units-per-pixel="0.08789062500000" order="3"/>
        <TileSet href="4" units-per-pixel="0.04394531250000" order="4"/>
        <TileSet href="5" units-per-pixel="0.02197265625000" order="5"/>
        <TileSet href="6" units-per-pixel="0.01098632812500" order="6"/>
        <TileSet href="7" units-per-pixel="0.00549316406250" order="7"/>
        <TileSet href="8" units-per-pixel="0.00274658203125" order="8"/>
      </TileSets>
    </TileMap>

openlayers.html:

var layer = new OpenLayers.Layer.TMS("Weather TMS","http://mysite-name.com/TMS/",
                    {
                         layername: "ne1",
                         type: "png",
                         isBaseLayer: true,
                         getURL:getURL
                    });
                   map.addLayers([layer]);

If there are no suggestions to fix this problem, could you point me to something similar to Tile Map Service. I am using gdal2tiles to break up my images into tiles and I have many images that I would like to serve through a tile service.


texts: How to change filtering by date for bimonthly GIMMS imagery to filtering by property 'year'?

texts:
In an earlier post at How to select features in a feature collection simultaneously row by row in Google Earth Engine?, I got help for creating a code to compare conflicts and bimonthly GIMMS imagery.
Now I received a different dataset, GIMMS annual maximum NDVI data, to use tor the time 1981-2015. The filtering from previously date ('system:time_start') which was added before, doesn't work now.
Also the NDVI (now called 'b1') band is now between -99999 to +9999.
//Dates for filtering the GIMMS Image Collection // With annual data no filtering nexessary, 1981-2015 used
var startDate = 1982 //adjust with time for the 35
var endDate = 2015
print(startDate);

//Loading and filtering (by date) the UCDP conflict data (CSV) 1989-2015
var conflicts = ee.FeatureCollection('users/welp1992/Conflicts_Deaths_Greater25')
  .filterMetadata('year', 'greater_than', 1989)//min 6 years, but as long as possible
  .filterMetadata('year', 'less_than', 2010) //including 2010, to 2015 = 6 years
print(conflicts.limit(1000),'conflicts');

//filter image GIMMS by date and band and add time 't' as a band
var timeSeries = ee.ImageCollection('users/welp1992/NDVI_Annual_Max_81_2015')
    .select(['b1'])
    .map(function (image) {
      return image.addBands(
        ee.Image(image.getNumber('system:time_start')) //'year'?
          .subtract(ee.Date(startDate).millis())
          .divide(1000*3600) // days (*365 = year)
          .int() //casts the input value to a signed 32-bit integer
          .rename('t')
      )
    })

print(timeSeries,'time');

This is the Error message I get:
ImageCollection (Error)
Error in map(ID=Max_05): Image.constant: Parameter 'value' is required.
Unfortunately I cannot share the annual max NDVI data (which I have uploaded to GEE). The script is basically the same as the earlier version in the link I shared above.
Do you have any advice?


texts: Updating ESRI mxd definition query variable via python?

texts:
we use definition queries to filter data to field users. We have up to 40-50 MXDs that are identical except for 1 variable in the definition query that allows the data for that MXD to be unique.
Example -  "DeploymentNumber='AA'" where 'AA' is the variable.
So, when we update MXDs, we would like to be able to use a python script to be able to update the variable - 'AA' component of our definition query in each layer without having to open each MXD, click on each layer and manually update (which in some cases is a 5 minute operation per MXD).


texts: Why is QGIS georeferencer much faster on Windows than Mac?

texts:
I just went from temporarily using a Windows 7 machine with 4GB of RAM and no video card as far as I know. I'm now using a mac with 16GB of ram, etc. Its a much more powerful machine.

However, QGis's georeferencer program runs about 8 times faster on the windows machine. Does anyone know why this is?

Gimp ran much faster too on the cheap windows machine.

Mac Specs:

2.5 GHz Intel Core i7
16 GB 1600 MHz DDR3
OSX Version 10.9.5
MacBook Pro Retina, 15-inch, Mid 2014
Solid State Hard Drive

Windows Specs:

4GB RAM
Windows 7
Data on External Hard Drive (non-flash)

I'm guessing its something related to the fact that videogames tend to run better on windows platforms, viz. more time has been put into the win32 software graphics engines to display real time graphics..


texts: How to calculate the center and radius of a circle based on points

texts:
I have a set of lat/long points. I need to calculate the center point and radius of a circle that will encompass all the points.

The co-ordinate system used is latitude/longitude as per google maps.

I previously used Think-Geo MapSuite where I could call a function, provide a collection of points and receive a bounding box. I would use the bounding box to instruct the map to display at a particular location and scale.

The tooling I'm using now is Xamarin Forms. This is a cross-platform tool that provides a high level wrapper across mobile platforms. For iOS this is Apple maps, Android is Google Maps and Windows Phone is Bing Maps. It is a new tool so the wrapper is very basic at present. I can only provide a center point and radius to instruct the map where to display.

I ultimately need C# code but pseudo code, or any other flavour would be a great help.


texts: QGIS composer add attribute table with repeating columns

texts:
I need some help with QGIS print composer.  I use the atlas composer and I try to use add attribute table from layout tab.

I have create my expression where I want to print in composer.

field expression (that I want):

case when "myfield" is not null   then concat(myfield,' ',myfield1,' & ',myfield2) else concat(myfield1,' ',myfield3)

I have defined filter from item properties where I want to print with atlas:

intersects(buffer($geometry),50),$atlasgeometry)

but that shows like this:

value1
value2
value3
value4
value5
value6
value7
value8
value9

because of to many values I can't print them all in my layout.  I need a result like this:

value 1   value2
value 3   value4
value 5   value6
value 7   value8
value9

If the values are like this then I can print all.  If I want to add a legend then have option where you can choose columns

That is what I need for my attributes but I can't find something like this.

Finally I try to create a new expression with my field and I try to add next feature but without success.

I tried this:

dbvaluebyid('test','f_field', @row_number ) ||'' ||dbvaluebyid('test','f_field', @row_number +1)

but from this expressions I don't take next feature.

and this:

case when "myfield" is not null   then concat(myfield,' ',myfield1,' & ',myfield2) else concat(myfield1,' ',myfield3)  end ||  '                                    ' || case when "myfield" is not null   then concat(myfield,' ',myfield1,' & ',myfield2) else concat(myfield1,' ',myfield3)

any idea?


texts: ESRI Leaflet not displaying simple Feature Layer

texts:
I'm having trouble getting my ArcServer feature layers to display at all in leaflet webmaps.  I'm just using a simple HTML tester, and I am able to display and mess around with ESRI-supplied feature layers, but when I link to my layers they simply do not display in the map.  I thought at first it might've been a projection issue as I was initially serving the data in local UTM, so I served up the data in web mercator to no avail.

  <meta charset=utf-8 />
  <title>Simple FeatureLayer</title>
  <meta name='viewport' content='initial-scale=1,maximum-scale=1,user-scalable=no' />

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.1/dist/leaflet.css"
    integrity="sha512-Rksm5RenBEKSKFjgI3a41vrjkw4EVPlJ3+OiI65vTjIdo9brlAacEuKOiQ5OFh7cOI1bkDwLqdLw3Zg0cRJAAQ=="
    crossorigin=""/>

I'm using this HTML tester to mess around with my code.  And this https://services.arcgis.com/rOo16HdIMeOBI4Mb/arcgis/rest/services/Heritage_Trees_Portland/FeatureServer/0 feature layer is the one supplied by ESRI that displays just fine when you switch it out for my service.

Any ideas?  Is there any specific formatting that data needs before it will display in leaflet?


texts: How to do better quality of texture on DEM in qgis2threejs?

texts:
I use as texture for DEM image with original resolution 4000x4000px. But when qgis2threejs creating 3D model in browser - texture looks like image with low quality.
In plugin menu there is line "resolution"

when i choose 200% or 400% quality are not changing.
How to do better quality of texture?


texts: Export 2D ArcScene into Adobe illustrator as vector data

texts:
I have an ArcScene (10.8) document with shapefile layers. I would like to export the 2d view into adobe illustrator for editing, but it appears as a raster. I suppose I could try to vectorize it in Ai, but that seems like a cumbersome workflow. This same layer, when exported from ArcMap, correctly shows up as an editable vector layer in Ai.
I see two similar questions from years ago that have not been answered:

Exporting 2D vector graphic from ArcScene to Illustrator
Exporting smoother contours using ArcScene

Does anyone know if this is possible?


texts: Transform point cloud from ENU (east, north, up) to LLA (lat, long, alt)

texts: I am having point cloud (xyz) in ENU format. I need a pipeline to convert all the points in the point cloud to geocentric LLA (lat, long, alt) coordinates. Is there any tool using which I  can achieve that?


texts: How to use a Value of an iterator row selection as a SQL request in ModelBuilder

texts:
I am trying to put as an input a XLS file with rows, each row represent the reference of an object.
I would like to iterate each row, and for each reference, collect all the data from my geodatabase.
Example: Reference = Object1, my ModelBuilder would search Object1 in my database and then use the data.
I found the way to get from the XLS to the value, but I don't know how to use the Value as an argument for my SQL request (expression).


texts: Don't add empty images to GEE ImageCollection

texts:
Given 2 dates I want to output all satellite images that 100% fill a given region (region is an imported shape file of a UK county boundary).
Initially I was just getting all images, but sometimes they would not fill the whole county and I could not figure out a method to make sure they fill the whole region with filters.
To solve this I decided to get weekly composite images instead. My reasoning was that if the satellite took a photo of half a county, it probably captured the other half after (I really didn't want to spend time meshing images together to get a full county images, so I thought the composite would be a quick fix).
Now my issue is as follows: when there are no images taken in a given week, my code adds an empty image to the image collection - I don't want it to do this, I only want images with actual content added to my collection - is there an easy way to do this?
My code is as follows:
function addNDVI(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']);
  return image.addBands(ndvi);
}

var clipArea = function(image){
  // table is my imported shapefile
  return image.clip(table)
}

function getWeeklySentinelComposite(date, product){
  var sentinelWeek = ee.ImageCollection(product)
  .filterBounds(table)
  .filterDate(date, date.advance(1, 'week'))
  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
  .map(clipArea)
  .map(addNDVI)
  .mean()
  .set('system:time_start', date.millis(), 'dateYMD', date.format('YYYY-MM-dd'));
  return sentinelWeek;
}

var startDate = "2018-01-01";
var endDate   = "2018-12-31";

var weekDifference = ee.Date(startDate).advance(1, 'week').millis().subtract(ee.Date(startDate).millis());
var listMap = ee.List.sequence(ee.Date(startDate).millis(), ee.Date(endDate).millis(), weekDifference);

var weeklySentinel = ee.ImageCollection.fromImages(listMap.map(
  function(dateMillis){
    var date = ee.Date(dateMillis);
    var image = getWeeklySentinelComposite(date, "COPERNICUS/S2");
    return image
  }
));

At this point weeklySentinel is an image collection of weekly image composites, some images are understandably empty. Is there a way to only add images with content to my collection, or do I have to try fortnightly / monthly composites and adjust my cloud cover filter to get a full dataset?


texts: How to find if polygon has a hole using field calculator in ArcGIS

texts:
Some 3rd party software cannot handle polygons with holes (donut polygons).
To find them I always use ArcView 3 field calculator, because:

however when I try to do this in ArcGIS, I get this:

So the question is how to make it right in ArcGIS.


texts: Does QGIS menu have any XY tool like in ArcGIS to add XY data? I have used the other method to capture XY coordinate. Can anyone help?


texts: Defining "Style file" of QGIS modeler process "Set style for raster layer"?

texts:
I want to add "Set style for raster layer" to QGIS processing modeler but I have no clue how I need to define the  "Style file" field. I have a layer style document and I assume, I need to define a path there. However adding the path to the directory works as little as adding only the layer style file-name. Drag'n'drop doesn't work either. Further it's not the gaps in the path, I tried it with a Path without spaces.

It's always the following error in the Log:


texts: How to print properties of a clicked feature to a panel in Earth Engine Apps

texts: I have an Earth Engine App where I display a ee.FeatureCollection. I want to be able to click on a feature and have its properties printed to a ui.Panel widget. How can I do this?


texts: Sorting by multiple fields at same time in QGIS?

texts:
Is it possible to sort multiple fields at the same time in QGIS?

For example, I'd like to sort a table by state name, and then within each state, sort by county name.

This would be the same as the Advanced Sort option available in ArcMap.


texts: Convert length in decimal degrees to feet

texts: I've got attributes for a length that the Near tool calculated, and it gave them to me in decimal degrees.  I need them in feet.  Is it possible to do this with the Calculate Field tool?  I tried reading the documentation in the ArcMap help, but it was really unhelpful.  I'm willing to do it in python or VBScript or whatever it takes.


texts: Editor Extension in ArcGIS Engine 10 Application

texts:
How can I find the Editor extension in an ArcGIS Engine application (i.e. no access to the IApplication object)?

There is an example method at http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//0001000004nn000000 ...

private void GetEditorExtension(IApplication application)
{
    //Find extension by CLSID or ProgId.
    UID extensionID = new UIDClass();
    extensionID.Value = "esriEditor.Editor";
        //{F8842F20-BB23-11D0-802B-0000F8037368}.
    IExtension editExtension = application.FindExtensionByCLSID(extensionID);
}

However, IApplication is only available for code that runs directly in ArcMap, not stand-alone apps.

If it helps, I'm writing this inside a custom tool, so I also have acces to the HookHelperClass.


texts: "ImportError: No module named gdalsrsinfo" What's the issue?

texts:
in the python console I can call import gdal, and gdal is recognized as a module; however, most gdal programs, such as gdalinfo and gdalsrsinfo are not recognized. I checked the GDAL source files to see if there was a problem with the PYTHONPATH. They were not in the QGIS application resources, so I moved them manually… this still did not solve the problem.

if i enter gdalsrsinfo /Filepath/File.prj the console returns the error:

  Traceback (most recent call last):   File "", line 1, in
     File
  "/Applications/QGIS.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 453, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level) ImportError: No module named gdalsrsinfo

However this is a function that should be included in GDAL by default. Why is it not recognized?

I am using the latest version from KyngChaos


texts: GRASS r.in.lidar tool generates incomplete rasters

texts:
I have been trying to import a huge LIDAR file (10.5) Gb to GRASS GIS using r.in.lidar. It does the task, but only for part of the image, as shown below:

Canopy Height Model (CHM) generated with GRASS:

CHM generated with FUSION:

The DEM, DSM and CHM generated in GRASS were built from the following commands, respectively:

r.in.lidar -o input=/my/route.las output=lidar_min method=min
r.in.lidar -o input=/my/route.las output=lidar_p99 method=percentile pth=99
r.mapcalc "CHM = lidar_p99@ig - lidar_min@ig"

All the three layers are incomplete, only some parts of the image are processed. I have tried different solutions: importing the files that I fused for creating the big one, importing files that had been created using FUSION (a file with only the ground points). I tried exporting the above rasters to geotiff and open them in QGIS, but the result is the same.

The only thing that worked is converting the LIDAR file to ASCII (only the xyz columns), and then importing it with r.in.xyz.

r.in.xyz input=/my/route.txt output=lidar_min separator=comma method=min

and it worked. This is the DEM produced using r.in.xyz from GRASS.

Then, the data are complete.

I think one possible reason behind r.in.lidar is not working as expected might be that GRASS only imports the pulses that are on the nadir or very near to it or only some passes of the flight. Does somebody know the reason behind this behavior?

 If it helps, this was the command line used in Fusion to generate the CHM that worked:
    CanopyModel /outlier:0,20 /ground:D:\DTM_default.dtm /ascii D:\CHM_default.dtm 5 m m 0 0 0 0 D:\laslist.txt

Sample

I uploaded a small sample of the data here (https://goo.gl/0S9p0l). It contains a file with only the xyz coordinates (.txt), other with the whole table (x,y,z, intensity, scan angle, etc.) (.las) and other two files that are DSMs created from the point files (.txt and .las).


texts: Converting raw encoded XML in KML into geojson properties

texts:
I have KML generated by a WMS (Web Map Service). I need this KML in geojson format and I'm doing this all on a linux box. My problem is that the description in the KML contains encoded HTML like so (some text changed):

<description>&lt;h4&gt;namehere&lt;/h4&gt;
&lt;ul class="textattributes"&gt;
  &lt;li&gt;&lt;strong&gt;&lt;span class="atr-name"&gt;id&lt;/span&gt;:&lt;/strong&gt; &lt;span class="atr-value"&gt;15L&lt;/span&gt;&lt;/li&gt;
  &lt;li&gt;&lt;strong&gt;&lt;span class="atr-name"&gt;issuetime&lt;/span&gt;:&lt;/strong&gt; &lt;span class="atr-value"&gt;Oct 14, 2016 8:39:56 AM&lt;/span&gt;&lt;/li&gt;
  &lt;li&gt;&lt;strong&gt;&lt;span class="atr-name"&gt;stormname&lt;/span&gt;:&lt;/strong&gt; &lt;span class="atr-value"&gt;NICOLE&lt;/span&gt;&lt;/li&gt;
  &lt;li&gt;&lt;strong&gt;&lt;span class="atr-name"&gt;validtime&lt;/span&gt;:&lt;/strong&gt; &lt;span class="atr-value"&gt;Oct 14, 2016 6:00:00 AM&lt;/span&gt;&lt;/li&gt;
  &lt;li&gt;&lt;strong&gt;&lt;span class="atr-name"&gt;windspeed&lt;/span&gt;:&lt;/strong&gt; &lt;span class="atr-value"&gt;64&lt;/span&gt;&lt;/li&gt;
&lt;/ul&gt;
</description>

Which looks like this after I run it through sed to fix the < and >:

<description>
namehere

id: 15L

issuetime: Oct 14, 2016 8:39:56 AM

stormname: NICOLE

validtime: Oct 14, 2016 6:00:00 AM

windspeed: 64

</description>

So my question is - is there an easy way to get those as properties into the geojson file that gets created when I run ogr2ogr (with ogr2ogr -f "GeoJSON" feed.json feed.kml) on it? My current plan is:

run it through sed to convert the < and the >
use xidel to grab each piece out of it with: "xidel feed.kml -e /kml/Document/Folder/Placemark/Description |sed '/^\s*$/d'" which will output multiple sections that look like this:

namehere
  id: 24W
  issuetime: Oct 14, 2016 8:39:56 AM
  stormname: NICOLE
  validtime: Oct 14, 2016 6:00:00 AM
  windspeed: 64

convert it to geojson with "togeojson feed.xml >feed.json"
change the description with jq
somehow loop and insert (using jq) the XML properties from the xidel output

This should work once I get all of these wired up but I can't help but think that there's an easier way. Is there some way to do all of this with a command line switch in ogr2ogr or via some other npm tool or something? My goal is to have the name/value pairs in each properties section, though honestly, if I could just get them somewhere accessible in the json I'd be happy.

UPDATE: I just realized I can use togeojson to convert the KML directly to json and that creates the json with the description containing valid HTML (with some stray \n strings). That doesn't require me to run it through sed. Maybe that would be a better path to pull out and change them into properties. Maybe. I may try that if nobody gives me any good ideas.


texts: Averaging large number of raster layers using ArcGIS for Desktop?

texts:
I have a large number of rasters (300+) I would like to average using Raster Calculator as is explained at Averaging set of rasters using Raster Calculator of ArcGIS Desktop? I know I can manually select each raster by double-clicking, hitting the '+' key, and then dividing by the total number of files.

Is there a less tedious way to get the average raster, without having to double-click and hit "+" hundreds of times? (or typing each raster's filename into the box, or something similarly tedious)


texts: Change size, color for nested styles consisting of different symbol layers created with geometry generator in QGIS

texts:
In QGIS 3.16, I have a nested linestyle that consists of different symbol layers. I am able to change the width, color and other settings manually for each symbol layer. However, I want to be able to change the setting for the whole composed symbol at once, as it is possible with some pre-defined styles.
In the example below (topo railway) is a pre-defined style that is shipped with QGIS and is for illustration purpose to show how the solution should work. The topo railway style consists of a marker line (no. 2) and a simple line (no. 3). If I change the size of the composed linestyle (no 1), the size of no. 2 and 3 changes proportionally.
How to achieve that if you start to compose an own style with different symbol layers?
Updated question: I have symbol layers created with geometry generator. See here for an example: https://gis.stackexchange.com/a/380267/88814
Especially these symbol layers are the problem - with other kind of symbol layers, setting the global size works.
Update 2: As a workaround, I can use a variable (global or per project) that can be used to define the size, width etc. of the different symbol layers. However, this is not very useful if the style thus created should be exported (using QGIS style manager) and used on another computer. One will have to define the variable separately on each system (or project).


texts: QGIS 3, OrfeoToolbox 7, Dynamic Library Error

texts:
I installed QGIS using KyingChaos’ OSX build and installed OTB following guidance on the GitHub page. When I run the kmeans classifier I get the following error:
Algorithm ‘KMeansClassification’ starting…
Input parameters:
{ ‘in’ : ‘/c1d_arc.tif’, ‘out’ : ‘/c1d_kmns.tif’, ‘nc’ : 5, ‘ts’ : 100, ‘maxit’ : 1000, ‘centroids.in’ : ‘’, ‘centroids.out’ : /centroids.out.file’, ‘sampler’ : ‘periodic’, ‘sampler.periodic.jitter’ : 0, ‘vm’ : None, ‘nodatalabel’ : 0, ‘cleanup’ : True, ‘rand’ : 0, ‘outputpixeltype’ : 5 }

dyld: Library not loaded: @rpath/libOTBCommandLine-7.0.1.dylib
Referenced from: /OTB-7.0.0-Darwin64/bin/otbApplicationLauncherCommandLine
Reason: image not found

GDAL was giving a similar error after I messed up the library installation, which I backtracked and fixed. So I’m thinking this is some similar sort of library version incompatibility. Also, in the readme it mentions being specific to Python 3.5, and I'm running Python 3.6.


texts:
When displaying two features, a polygon a (see geoJSON) and its minimum bounding rectangle b (see geoJSON), I get the following result when viewing on geojson.io:

The minimum bounding rectangle around the polygon looks more like a parallelogram than a rectangle to me.

However, when I upload these same features to OpenLayers, the distortion is gone:

Is this a projection issue?


texts:
I'm trying to follow this tutoriel to do a concave hull script :
Drawing Boundaries In Python
Of course, I drop all steps for plot results. But at the and of my script I have an error.
This is last lines :
edges = set()
edge_points = []

# loop over triangles:
# ia, ib, ic = indices of corner points of the triangle

for ia, ib, ic in tri.vertices:
    pa = coords[ia]
    pb = coords[ib]
    pc = coords[ic]
    # Lengths of sides of triangle
    a = math.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
    b = math.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
    c = math.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
    # Semiperimeter of triangle
    s = (a + b + c)/2.0
    # Area of triangle by Heron's formula
    area = math.sqrt(s*(s-a)*(s-b)*(s-c))
    circum_r = a*b*c/(4.0*area)
    # Here's the radius filter.
    #print circum_r
    if circum_r < 1.0/alpha:
        add_edge(edges, edge_points, coords, ia, ib)
        add_edge(edges, edge_points, coords, ib, ic)
        add_edge(edges, edge_points, coords, ic, ia)

m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points

And it returns :
Assertion failed: (!static_cast<bool>("should never be reached")),
function itemsTree, file AbstractSTRtree.cpp, line 373.
Abort trap: 6

What is wrong here?


texts: Setting color of Simple Marker Layer from field using python 3

texts:
I'm trying to set up the colour of a QgsSimplerMarkerSymbolerLayer.

Here is my code :

sites=iface.activeLayer()
single_symbol_renderer=sites.renderer()
symbol=single_symbol_renderer.symbol()
symbol.setSize(4)
prop = QgsProperty()
prop.setField("Couleur")
pc = QgsPropertyCollection('Color')
pc.setProperty(4, prop)
marker_symb = QgsSimpleMarkerSymbolLayer()
marker_symb.setDataDefinedProperties(pc)
sites=QgsSimpleMarkerSymbolLayer(marker_symb)

But it still not working, any help ?


texts: How to dissolve polygons and aggregate their data in QGIS?

texts:
I have a shapefile (one layer) with my country but I have a problem.

It is divided into 4 entries (prefectures), and I want it as one, not 4.

Also, there is a table with some data, and I want to have one row with the sum of the data. For example:

Entry  1: population 10.000,
entry  2: population 20.000,
entry  3: population 30.000,
entry 4: population 40.000.

Ultimately I want something similar to  this:

main piece: population 100.000.

(I don't care if this will be done by hand, its not much work!).

I'm a new QGIS user and i don't know much. Any help would be really appreciated greatly!


texts: Google Earth Engine: Masking land out of Sentinel-3 imagery

texts:
I'm trying to isolate ocean cover from OLCI imagery using a bitmask from the quality_flags band. But it's masking all the pixels, and I notice when I mask the ocean and leave the land, the land is very patchy compared to non-masked imagery.

Here's a link: https://code.earthengine.google.com/8b979aedb7f91352ccd792c63af8ca39

Here's the code I'm using:

//Land mask
function maskLand(image) {
  // Bit 31 is land/sea cover
  var LandBitMask = ee.Number(2).pow(31).int();
  // Get the pixel QA band
  var qf = image.select('quality_flags');
  // The flags should be set to one, indicating water cover.
  var mask = qf.bitwiseAnd(LandBitMask).eq(1);
  return image.updateMask(mask)
}

//Create image collection for area and date range
var lis_collection = ee.ImageCollection(olci)
  .filterBounds(lis)
  .filterDate('2018-05-01', '2018-10-01');

//Mask the image collection to water
var lis_mask = lis_collection.map(maskLand);

//Mosaic, clip, and display the area
var lis_area = lis_mask.mosaic();
var lis_area = lis_area.clip(lis);
print(lis_area);
Map.addLayer(lis_area, oa5vis, "Long Island Sound radiance");

Is there a better way to do this?

Perhaps some kind of mask I can run on the clipped mosaic instead of a function to run on the image collection.


texts: ADC Maps Grid to Lat/Long Coords Conversion

texts: While the ADC site apparently has a paid conversion function, it's interactive, and I'm looking for some way to do so on a batch basis, possibly with an API.  Any hope?


texts: proj4js, proj4.net, proj4 give different result to cs2cs

texts:
I'm trying to convert between two arbitrary projected coordinate reference systems, but every proj.4 library variant that I've tried has given a results different from QGIS's cs2cs application (which matches EPSG.io's results).

So, for transforming the point (572084.7402, 4542311.793) metres, the expected result is (980191.1775, 1086873.446), but the libraries return (980156.3975, 1087102.081).
The jsfiddle is here: http://jsfiddle.net/dougy83/nfr62bj1/2/ , which includes the proj4 definitions for the reference systems.

I've used are the JS, C and .NET versions of the proj.4 library, and they all give the wrong value.

I've used the QGIS cs2cs, and EPSG.io and mygeodata.cloud/cs2cs/ websites, which all give what I assume is the correct answer. The cs2cs code appears to just call the pj_transform function, which is all I'm doing for the C version library.

Could someone please provide some help with this? I keep seeing references to the datum shift / grid shift files; is this related to this problem?

Proj4 definitions are duplicated here:

source: +proj=lcc +lat_1=30.28333333333333 +lat_2=28.38333333333333 +lat_0=27.83333333333333 +lon_0=-99 +x_0=600000 +y_0=3999999.9998984 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs
dest: +proj=lcc +lat_1=27.41666666666667 +lat_2=34.91666666666666 +lat_0=31.16666666666667 +lon_0=-100 +x_0=914400 +y_0=914400 +datum=NAD27 +units=ft +no_defs


texts: How do I convert a string into a column number in QGIS Expression Builder?

texts:
I want to highlight parts of a map based on the page number in an atlas. Ideally I'd use the expression builder in the layer style menu to set the Column to equal the atlas page number.

My expression is as follows

@atlas_pagename  || '_1'

The page number is "1007" and the column is "1007_1" hence the need to add the _1 at the end. The next column is 1008_1 for example.

The expression successfully outputs 1007_1, but the expression builder doesn't recognise this as the name of a column.

Is there any way to make the expression builder not just output "1007_1" as a string, but to select the column "1007_1"?

I am trying to do this automatically with an expression so that when the atlas moves onto page 1008 the layer will move onto column 1008_1 and so on.


texts:
I have some census tract polygons with language spoken at home attribute data behind them. I created them in ArcMap and produced a dot density map. However, after importing to AGOL it's only the polygons. I am aiming to create a swipe map comparing two maps visualizing the growth of language outwards from the city to suburbs over 20 years.

Here's a before and after image of the tracts in ArcMap, then AGOL:


texts: MMQGIS Plugin dropdown/icon not appear (QGIS 3.4)

texts: I am on the latest LTR for QGIS (3.4) and have installed the MMQGIS plugin. In the older version of QGIS (2.18) this caused a dropdown menu to appear on the top menu bar. However, there is no dropdown menu or icon in QGIS 3.4 in order to access the plugin features. I cannot figure out how to access the MMQGIS plugin short of going into the python repository for it. Any suggestions on how to manually add a link to the plugin or address the issue?


texts: how to define the corrected feature extent after using conversion tools in arcgis 10.2?

texts:
I have a map in dgn format. After using conversion tools in arcgis to extract river layer, the feature extent as followed

and the project is unknown. I can define the project using define project in arctool box, however, the feature extent is weird, I want to change to the corrected one.
the corrected extent should be

How do I redefine the feature extent of that map?


texts: OSGB feaures not aligning with OpenStreetMap basemap even though Data Frame and Features are all in OSGB?

texts:
I have a task to report features in OSGB (British Nat Grid) over an OpenStreetMap basemap.

I have done this successfully in past with line up of features on basemap, but this time I cannot get the features correctly aligned.

Comparing the two (old test and new task) map files, I cannot see any difference in settings but one works and one does not. All features (vectors) are in OSGB and the dataframe for both is in OSGB, but something else must be causing the offset.

Anything subtle I might have missed in this new task.

I did spend a lot of time in the previous test exploring this but can't remember anything that seems missing now.

I have tried for the new set, creating new maps and adding alternately the base map and features first and in both cases making sure data frame is OSGB, but no difference to the result.


texts: Geopandas sjoin operations fail to capture only polygons within other as it excludes points on the boundary

texts:
This is the code I'm using:
result = geopandas.sjoin(tracts, city, how='inner', op=X)
(where X is the operation I'm changing).
From the documentation these are some of the options for op and the results I've got:

covered_by: "Returns True if every point of object is a point on the interior or boundary of other."

intersects: "Returns True if the boundary or interior of the object intersect in any way with those of the other."

overlaps: "Returns True if the geometries have more than one but not all points in common"

within: "Returns True if the object’s boundary and interior intersect only with the interior of the other"

I need to get all of the tracts within the city (and only within the city, including the ones that touch the edges). From the documentation, I thought that covered_by would give me this, but apparently not. What op should I be using?


texts: Instruction to import OSM file using PostGIS in windows

texts: I want to import OSM data into a PostGIS database running under windows. After a web search I learned that I should use osmosis. But most of the explanations are for Linux. But I have to use windows. Is there a detailed description how to do that with Windows?


texts: ArcGIS Pro rotate point symbol smaller than one degree

texts:
Is it possible to rotate a point symbol by an angle smaller than one degree?

When i try to enter 0,5° in the Angle field of the Symbol pane it gets rounded to 1°. 0°30′ doesn't work either.

My workaround solution would be to rotate the underlying .svg file, but that's not a nice solution.

Maybe it is feasible in Python?


texts: Change simple symbol color by expression

texts:
Using QGIS, I want to display points on a basemap, the color of the point (symbol) I need to change from a 'rule based' expression.

I have loaded a basemap, imported my data, a .csv file containing amongst others a lat column, a lon column and a 'level' column, containing values from -30 to -140.

I have selected 'rule based' and in the 'expression' I have "level" >-108.
This puts me a colored dot on the map at the correct point when the value in the 'level field is greater than -108, but what I want is for the dot to be a different color if 'level is less than -108, (there are some other fields that I will deal with later, hence needing a rules based expression).

So, basically, my question is how do I make the symbol a different color depending on the value when using expressions? I have tried color_rgb and symbol_color, neither worked.


texts: Running custom model via ArcPy with PyCharm fails with exit code -1073741819 (0xC0000005)?

texts:
I'm trying to geocode a few address using a custom model created. When running the model via model builder all runs perfectly.

Also when running it via the embedded interactive python console in ArcCatalog using the following code it also runs smooth.

>>> arcpy.ImportToolbox("C:\\Users\\jvidin\\AppData\\Roaming\\ESRI\\Desktop10.1\\ArcToolbox\\My Toolboxes\\usps_geocoder.tbx","geocoder")
<module 'geocoder' (built-in)>
>>> arcpy.geocoder_geocoder(
... 'gc_01',
... 'J:/Postal/Postal/99_Geocode/01/composites/composite_1',
... 'J:/Postal/Postal/99_Geocode/01/composites/composite_2',
... 'J:/Postal/Postal/99_Geocode/01/gc.gdb/gc_01'
... )
...
<Result ''>
>>>

The problem occurs when trying to run it via IDE using the following code:

import arcpy

arcpy.ImportToolbox("C:\\Users\\jvidin\\AppData\\Roaming\\ESRI\\Desktop10.1\\ArcToolbox\\My Toolboxes\\usps_geocoder.tbx",
                    "geocoder")
try:
    arcpy.geocoder_geocoder('gc_01', 'J:/Postal/Postal/99_Geocode/01/composites/composite_1',
                        'J:/Postal/Postal/99_Geocode/01/composites/composite_2',
                        'J:/Postal/Postal/99_Geocode/01/gc.gdb/gc_01')
except Exception as e:
    print e.message

Generating the following error:

  *C:\Users\jvidin\2.7_32_arcpy_pyenv\Scripts\python.exe C:/Users/jvidin/PyCharm_projects/2.7_ArcPy_Rebuild_locators/run_geocoding_model_standalone_stackexchange.py

  Process finished with exit code -1073741819 (0xC0000005)*

Model created in arc 10.1 and running with arcpy 32 bits 10.1 python install in a virtual env.

## new edit 12/15/2014 . continuation ...

So i changed a bit my model and my script, this time the model has one and only one exposed parameter and its an inline variable %var% used for outputs paths creation .
Script below.

  it creates the output .gdb based on argument sbr and sets env.workspace to it. (works)
  than when it runs arcpy.geocoder2_geocoder(sbr), in egg. sbr = '01' it fails with:
  Process finished with exit code -1073741819 (0xC0000005) not of much help.

def geocoder(sbr):
    try:
        fgdb = str(sbr) + '.gdb'
        path = "J:/Postal/Postal/99_Geocode/_results/"
        env.workspace = path+fgdb
        print env.workspace
        arcpy.CreateFileGDB_management(path, fgdb)
        arcpy.geocoder2_geocoder(sbr)
    except Exception as e:
        print e.message

really have not clue why it continues to fail, all the outputs of geocode2_geocode() model are setted to same env.workspace created.

The geocode2_geocode() model -> contains 2 composite locators that run one after other.


texts: QGIS TIN interpolation - error in pixel size of output

texts:
When using the "TIN Interpolation" tool in QGIS 3.10, the output raster file doesn't keep the pixel size I specify. For example if I specify a pixel size of 1.0m, the output file actually has a pixel size around 0.99m. The pixels also aren't perfectly square sometimes. Has anyone else experienced this problem or discovered a solution?

I have pasted below example screenshots of two parameter values, and the properties of the output files I get.

Overall for this example I want to get an output raster with 500 rows and 100 columns, pixel size 1.0m. It's impossible to enter this in the parameters because when I change one value, another changes. e.g. when I enter a pixel size of 1.0m, the 'Rows' parameter updates to 500 and the 'Columns' parameter updates to 101. If I change the 'Columns' parameter to 100, the 'Rows' parameter automatically changes to 495 and the pixel size changes to 1.0101.


texts: Why does my Google Earth Engine code run on a polygon, but not an image if they are of the same area

texts:
I have an Image (TIFF) that I imported and would like to run this script over to calculate average NDVI over that region, but it does not work giving me an error message stating "ImageCollection (Error)
Feature, argument 'geometry': Invalid type. Expected type: Geometry. Actual type: Image<[b1]>."
However, if I manually create a polygon geometry over the exact same region and have that as my input to filterbounds over it works perfectly.
Why does it accept the "geometry", but not the "image" if it is of the same region?  In the example code below, I have the "image" as the filterbounds and it does not run as such.  I realize in the error message that it says it is expecting a "geometry", but I don't know how to get it to expect an "image".
code example:
//This approach starts with a collection of L8 TOA images and filters them by cloud cover

// --------------Set Location and Filter by Date and Cloud-------------------//
var MAX_LIST_SIZE = 100;

var L8_TOA = ee.ImageCollection("LANDSAT/LC08/C01/T1_RT_TOA")
    .filterDate('2013-12-25', '2016-12-25') //For whole collection, comment this out
    .filter(ee.Filter.eq('WRS_PATH', 41))
    .filter(ee.Filter.eq('WRS_ROW', 36)); //This eliminates multiple WRS results
var count = L8_TOA.size();
print('Collection size after date filtering, if that was set: ', count);

var boundary = ee.Geometry.Rectangle(-122.45, 37.74, -122.4, 37.8);

        // **** This is not working correctly bc it only allows for "geometry"
var c = L8_TOA.filterBounds(image);          // Geometry, Image, or Boundary ***

var withCloudiness = c.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: image,                         // Geometry, Image, or Boundary ***
    scale: 30,  });
  return image.set(cloudiness);});

// Here is where the computed cloud score threshold is set!
var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 5));
 print('L8 TOA Filtered Collection, by date and cloud',filteredCollection);


texts: How to create tiles with manage tile cache tool in ArcGIS Desktop from PNGs?

texts:
I have a directory of PNGs without any geographic information other than quadkeys which I am converting manually via Python to get the bounding box coordiantes.  For testing purposes I would like to create a single tile in the location that is associated with it.  The steps that I have taken are, create mosaic dataset, apply the extent (left=-10350185.126039146, bottom=4315940.3650941895, right=-10348962.133586584, top=4317163.357546754) to the mosaic dataset, tile add individual raster with above coordinate information, then manage tile scheme.

The issues that I am seeing is that the extent properties in the layer and mosaic dataset are auto updating to 0, 255, etc.

What steps need to be taken to visualize the PNG as a managed tile cache with a mosaic data set and my specified extent?

Where should my processing extents be set so that the PNG that I am loading has the correct extent(around Louisiana) it is currently in Australia?  I am using web mercator auxillary sphere for my projection.


texts: Selecting features with NULL attributes in QGIS

texts:
How do I identify which features have NULL attributes in the attribute table of a shapefile?

When i try [field]= 'NULL' it selects the features with value = 0 (see picture)


texts: How do you insert values into custom fields when adding to a Mosaic Dataset table?

texts:
I have a mosaic dataset (file geodatabase) with an additional field added to it (id that is integer) and I'm trying to accomplish in ArcPy what I can do manually--which is:

Add a new raster (single tif file) to the mosaic dataset in the file geodatabase.
Set the additional field to a value.

I've successfully added a raster using  arcpy.AddRastersToMosaicDataset_management, but I don't see how I can set that additional field.

If AddRastersToMosaicDataset_management returned objectids, I could manually set it via an update cursor but I don't see how to get that value.

Is it possible from ArcPy to update a custom field on a Mosaic Dataset?


texts: How to speed up query performance with PostGIS data served up in GeoServer and OpenLayers?

texts:
I have published a table with @ 6 million polygons to geoserver, and have built a tool using the OpenLayers filter class to query the data. The performance is abysmally slow, often times not ever completing the query. I have indexed the columns I am filtering on and have also partitioned the table. The biggest issue is with the LIKE queries, and most recently I have applied the gin index as referenced in this post:

 https://blog.ianturton.com/postgis/2018/04/04/speeding-up-likes.html

But the fetch request with the query just spins forever when a LIKE filter is applied. On the desktop, the same queries run and fairly fast. When applying equivocal filters the query does execute, though still slower than in the desktop. What could be causing the slow down?

Here is the code I am using to generate the request, the filter_val variable is set via a form:

var featureRequest = new ol.format.WFS().writeGetFeature({
        srsName: 'EPSG:3857',
        featureNS: 'http://geoserver.sf.net',
        featurePrefix: 'sde',
        featureTypes: ['parcels_slim_part_view'],
        outputFormat: 'application/json',
        filter: filter_val
        // count: 100
      });

    controller = new AbortController();
    signal = controller.signal;

    var throbberCounter = 0;
    var increaseThrobberCounter = function() {
        document.getElementById('throbber').style.display = 'block';
        ++throbberCounter;
    };
    var decreaseThrobberCounter = function() {
        --throbberCounter;
        if (throbberCounter <= 0) {
            throbberCounter = 0;
            document.getElementById('throbber').style.display = 'none';
        }
    };

      // then post the request and add the received features to a layer

      increaseThrobberCounter();
      this.throbberVisible = true;
      fetch('http://localhost:8080/geoserver/wfs', {
        method: 'POST',
        signal: signal,
        body: new XMLSerializer().serializeToString(featureRequest)
      }).then(function(response) {
          console.log(response);
          console.log(response.url);
        return response.json();
      }).then(function(json) {
        var features = new ol.format.GeoJSON().readFeatures(json);
        if (features.length == 0){
            decreaseThrobberCounter();
            this.throbberVisible = false;
            alert('no records found, try again');
        }
        else{
        filterSource.addFeatures(features);
        map.getView().fit(filterSource.getExtent());
        decreaseThrobberCounter();
        this.throbberVisible = false;
        map.addLayer(filter_lyr);
        filter_lyr.setVisible(true);
        }
      });


texts: Grouping pixels into clusters based on pixel value in ArcGIS Desktop?

texts:
I have a raster image with different pixel values, each pixel values represent number of persons living in that pixel. I want all the pixels with population which are nearby to be clustered into one group.

How can i do this?

The image which should be clustered looks like this:


texts: Loading style (*.qml file) to layer via PyQGIS?

texts:
Is there any way to load a *.qml file to layer and set a style to it in python (pyqgis)?

Can I 'get into' this file in console? I want to review it.


texts: Exporting geodata from QGIS to Esri File Geodatabase?

texts:
For a project I need to create a set of vector data that I think I could do in QGIS.
However the customer demands that the final data set needs to be delivered in a Esri File Geodatabase.

Is there a way to export geodata from QGIS into an Esri File Geodatabase?

(I have only little experience with QGIS and not found a solution to this problem yet)


texts: Add cloud layers to vts-browser-js 3D map

texts:
I want to add 3D clouds to an existing 3D map based on the vts-browser-js library to simulate weather conditions.
I've added 3D 'obj' models to the map using the API.
Also vts-browser-js has a built-in skydome to show the air around the earth.
I wonder is there a way to add 3d clouds to the map?
Thanks in advance


texts: SQL Server UNion of geography feature returns FULLGLOBE result?

texts:
I have some trouble with sql server 2012.

I store coordinates of french cities in my database as geography object (with ESRI 4326) and i want to make the UNION of some cities.

I have in result (with the method AsTEXT() ) :

If I ask for 100 cities : FULLGLOBE
If I ask for 1 cities : I have the good coordinates

There is a limit of how many cities i can union ?

(I do my sql request in c# with ASP .NET MVC 4 :

DbGeography result = System.Data.Spatial.DbGeography.PointFromText("POINT EMPTY", 4326);

//PolygonFromText("POLYGON((0 0,0 1,1 0, 0 0))", 4326);;
string resultString = "";

try
{
  var commune = db.Communes
    .Where(com => com.IDCommune.StartsWith("073"))
    .Select(com => com);

  foreach(var com in commune)
  {
    if (com.coordinates == null)
    {

    }
    else
    {
      result = result.Union(com.coordinates);
    }
  }

Thank you !

EDIT / I solved my problem : i put my polygon in geometry type and not in geography type ans it works =)


texts: How to convert from HDF5 to GeoTIFF with Python

texts:
I'm coming from Reading metadata in HDF5 file.
I have a HDF5 file with some radar information. I'd like to read this data and geospatial metadata contained in the file to produce a GeoTIFF. Using the accepted answer, I can read the metadata, but I'm not sure how to convert it to something that rasterio can understand to produce a valid CRS and transform that allow to georeference the field.
Right now, the metadata I can get from the HDF5 are these two dicts:
{'projection_indication': b'Y',
 'projection_name': b'STEREOGRAPHIC',
 'projection_proj4_params': b'+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0'}

{'geo_column_offset': array([0.], dtype=float32),
 'geo_dim_pixel': b'KM,KM',
 'geo_number_columns': array([700], dtype=int32),
 'geo_number_rows': array([765], dtype=int32),
 'geo_par_pixel': b'X,Y',
 'geo_pixel_def': b'LU',
 'geo_pixel_size_x': array([1.0000035], dtype=float32),
 'geo_pixel_size_y': array([-1.0000048], dtype=float32),
 'geo_product_corners': array([ 0.      , 49.362064,  0.      , 55.973602, 10.856453, 55.388973,
         9.0093  , 48.8953  ], dtype=float32),
 'geo_row_offset': array([3649.982], dtype=float32)}

From the former I can easily obtain the CRS, but I don't fully understand the meaning of the parameters in the latter, although I'm pretty sure they implicitly define the transform.
How can this be used to produce a regular CRS and tranform understandable by rasterio?


texts: QGIS 3 processing script won't show gui when executed from editor

texts:
Different from QGIS 2, when I execute a processing script from the Processing Script Editor by pressing the Play-button in QGIS 3 it's GUI is not created. Seems to me that the initAlgorithm() method is not called?

Everything's fine with the GUI when double clicking the same script it in the processing toolbox:

Am I missing something? Issue?


texts:
I downloaded Cesium Java Script library for creating 3D Globes from "http://cesium.agi.com/".
Also I found that, they are using Bing satellite map as default. But I need to display layer from My Geoserver. I am very new to this Java Script library.

Can you please tell me , How to add a WMS Layer from Geoserver on the Cesium globe?
Is it possible?

Thanks


texts: Batch geocode city/state to lat/long using free or open source tool

texts:
I am looking for a way to batch geocode a dataset containing cites/states to lat/longs. The lat/long for the centroid of the city is fine. My data does not contain streets or street numbers. Only city/state combinations, e.g. Brooklyn, NY.

I have tried various geocoders, such as the Census Geocoder and the Google Sheets "Awesome Table" plugin "Geocode." Both of these geocoders require a "full address" and will not geocode city/state alone without a street/street number.


texts: How to calculate next point given two gps coordinates and speed?

texts:
If an entity is traveling in a straight line and at a known speed of 1km/min, how can I calculate the next point at time t given two gps locations? The altitude does not matter. The gps format is in EPSG:4326.

Can someone provide an example in python? It's not as simple as finding the slope of a line, like how we learned in algebra class?


texts: Importing GeoJSON to PostGIS with ogr2ogr and keeping geometry?

texts:
I have a .geojson file with line types that I'm trying to import into PostGIS with ogr2ogr. Following answers to similar questions on here, I tried the following:

ogr2ogr -f PostgreSQL PG:"dbname=gis user=alex" roads.geojson -nln roads

Using pgAdmin I can see that the roads table was created and the data imported, but when I go to add the table in QGIS in shows up as NoGeometry. I've tried adding -nlt with LINESTRING, MULTILINESTRING, and GEOMETRY, but no luck.

I'm new to PostGIS and ogr2ogr--can someone please shed some light on what I might be doing wrong?


texts: How to retain all attribute data when clip two shapefile?

texts:
I have two shapefiles, one is soil map and another one is paddock map.
I want to label every paddocks with soil type, but I don't know how.
I tried to use intersect , but it created more polygon like this.
Is there any options I can try?


texts: Geoserver Dynamic Raster Masking Based on Request Parameter

texts:
Is it possible for Geoserver to dynamically mask out values above or below a certain threshold in a WMS or WCS Geo-tiff request?

For instance I make a request with some threshold parameter like 49 and all values above that get masked out, set to 0, or -9999 (the no data value) in the response.

EDIT -- In response to BRADHARD's comment

With his suggestion I use the following request and SLD as a test:

myserver/geoserver/wms?service=WMS&request=GetMap&bbox=-125.020833333333,24.0625,-66.479166666662,49.937500000002&srs=EPSG:4269&env=opacity3k:0&layers=gdd:agdd&styles=gdd:agdd_web,&format_options=layout:provisional&width=1800&height=800&time=2016-06-17&format=image/png

-

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">
<NamedLayer>
<Name>agdd</Name>
<UserStyle>
<Name>agdd</Name>
<Title>AGDD</Title>
<FeatureTypeStyle>
<Rule>
<RasterSymbolizer>
<Opacity>1.0</Opacity>
<ColorMap>
  <ColorMapEntry color="#FFFFFF" quantity="0" label=""/>
  <ColorMapEntry color="#5e4fa2" quantity="1000" label="1000"/>
  <ColorMapEntry color="#3288bd" quantity="2000" label="2000"/>
  <ColorMapEntry color="#66c2a5" quantity="3000" opacity="${env('opacity3k',1)}" label="3000"/>
  <ColorMapEntry color="#abdda4" quantity="4000" label="4000"/>
  <ColorMapEntry color="#e6f598" quantity="5000" label="5000"/>
  <ColorMapEntry color="#ffffbf" quantity="6000" label="6000"/>
  <ColorMapEntry color="#fee08b" quantity="7000" label="7000"/>
  <ColorMapEntry color="#fdae61" quantity="8000" label="8000"/>
  <ColorMapEntry color="#f46d43" quantity="9000" label="9000"/>
  <ColorMapEntry color="#d53e4f" quantity="10000" label="10000"/>
  <ColorMapEntry color="#9e0142" quantity="11000" label="11000"/>
</ColorMap>
</RasterSymbolizer>
</Rule>
</FeatureTypeStyle>
</UserStyle>
</NamedLayer>
</StyledLayerDescriptor>


texts: How to equalise colour, brightness and contrast across several scanned paper maps

texts:
I am a vectors guy mostly, but for a current project I need to work with a bunch of old scanned paper maps (ww2 bomb damage maps for London, if your interested!)

We scanned and georeferenced the maps, and now want to create a tiled composite layer for serving on a website. We will cut out the borders obviously, thats not the issue here.

Problem is, that there are unsightly visual colour and brightness differences between map sheets. I am a bit at a loss as to how to equalise them so to give a coherent visual look. I read up on histogram equalising, but my current toolbox (Manifold GIS, GDAL, GeoServer) doesnt seem to have the required functionality.

Example of 4 scans already georeferenced:


texts: From this morning if I open a qgis project all the shapefiles are not editable anymore. I've tried different projects in both 2.14 and 2.16 qgis versions. Has anyone encountered the same problem?


texts: Is it possible to integrate CA SiteMinder with ArcGIS Server 10?

texts:
I am researching whether or not it is possible and/or desirable to integrate CA SiteMinder with ArcGIS Server 10 in order to restrict access to the various web service endpoints based for authenticated users.

I have googled for whitepapers, solutions, and documentation but every search is coming up blank.

Has anyone had experience with a CA SiteMinder / ArcGIS Server integration?
If so, is there any documentation on the integration?


texts:
Are there any tools that can take a DTM, in asc or geotiff format and generate a mesh from this that can used in meshlab or similar software.

I have asc DTM file that I need to convert to a ply mesh file. Is this possible?

In response to the comment. Then these DTM asc files are not the ones that meshlab understands.

The files starts with:

ncols        6250
nrows        6250
xllcorner    630000.000000000000
yllcorner    6070000.000000000000
cellsize     1.600000000000
NODATA_value  -9999


texts: How to visually sort the order of Points with QGIS, GDAL, ...?

texts:
The following Points layer is the basis of a rendering in a Android application.

Android allows to draw a shape from the provided points. I cannot use LineStrings.

However, the order of points exported from QGIS seems to be odd.

Is there any way to visually sort the order of points in QGIS or using the GDAL tools? The order of feature ids in the final export format (such as CSV) needs to be kind of a roundtrip as visualized in the following picture. If there are multiple paths or a "hole" it is perfectly okay to have multiple point layers.

Export experiments

I tried different steps to get hold of the latitude/longitude coordinates:

QGIS: Export points layer as CSV, add layer create options: GEOMETRY=AS_YX
Convert into Android specific format JSON format using regular expressions.

This has the disadvantage that the feature id is still included. The order of points however is not as desired.

ogr2ogr -f CSV output.csv input.shp -lco GEOMETRY=AS_YX
Convert into Android specific format JSON format using regular expressions.

Here the feature id is already omitted. The order of points however is not as desired.

QGIS: Export points layer as GML
ogr2ogr -f "CSV" -select "" output.csv input.gml -lco "GEOMETRY=AS_YX"
Convert into Android specific format JSON format using regular expressions.

The order of points however is not as desired. The order of points is different to what the former exports/conversions output.

The original layer is also available as a LineString layer.


texts: Summary Statistics after Spatial Joining in Modelbuilder

texts:
I am currently building a model to help speed up part of the generation of a map I make occasionally for the Dept of Public Works I work with. I have the majority of the map already made and saved in an .mxd. What I'd like to do is be able to open the model, run it and it adds the extra two new layers built from other existing layers and then have the completed map.

Manually what is being done is this: I have the point layer Water_Breaks_1997_to_Current, which I run a Select by Attribute SQL search to remove approx 100 results that are in the layer but no longer needing to be mapped. This result is then saved as a new feature class Water_Breaks_on_Active_Line_DATEHERE. Then right clicking my other layer representing water mains know as wMain (line geometry), I execute a spatial join where the layer wMain to the newly made Water_Breaks_on_Active_Line_DATEHERE. A points to line with one to many cardinality. Then summary statistic it by having each line intersecting a point(s) get a SUM count of said points, so that I can symbolize the new feature class to show the count of which water mains have the most breaks occurring on it in the city.

The problem that is occurring in the model is the spatial join. When opened, it does not offer spatial statistics the same way as right clicking the feature class itself. But through searching older SE questions and geonet, I came up with you have to add a new field, right click for merge rule and I select the 'Sum' ability.

At this point, I've validated the model, and run it in Edit mode multiple times. Each time it does everything I wish...except the summary statistic of the amount of points on each line. It shoots out a field populated full of 0 each time.

What step or option am I missing that is stopping this from being counted properly?

Using: ArcGIS 10.2.2


texts: Arcpy CalculateField using formatting syntax

texts:
I have 2 feature classes (FC) that I am merging together and then performing some field calculations on.

The 1st FC has 100k+ records in it, and the second FC is a blank FC that I'm using as a template. I need to merge the two FC's together and then field calc the values from FC1 into FC2.
The problem with FC2 is that the fields are required to often be of a different type. For example, in FC1 field "PID" is a double and in FC2 "PIN" is a string. I need to calculate the values of "PID" into "PIN" with a prefix of "139-". This pattern is repeated several times with different fields.

My simple code example:

data = r"C:\Temp\Test.gdb\PAR_PROPINFO_SCO_VW2"
template = r"C:\Temp\Test.gdb\Template"
out = r"C:\Temp\Test.gdb\Scott_CO"
arcpy.Merge_management([data,template],out)
qry= '"139-"' + "+'PID'"
arcpy.CalculateField_management(out,"PIN",qry)

In this particular instance, qry threw a syntax error.

If I just use:

arcpy.CalculateField_management(out,"PIN","PID")

I get a type mismatch error.

Since I'm on Arc 10.0, I was trying to avoid a cursor with the amount of records involved.

This process will run daily and the # of records will only increase. I was thinking I could just get away with a brute force merge and calc method. However, I run into type mismatches. I thought field mappings, but that is beyond my current level.


texts: Android offline app to visualize GeoTIFF maps (or any well-known raster format)

texts:
I like to go to mountain, and got lots of GeoTIFF maps, which I can convert to other formats with GlobalMapper (JPG+JGW, ECW, ...).

However, I haven't found an Android app to:
- show these raster images
- show my position, at the same time

I tested QGIS Mobile, but it did not show my position. Any other alternative?


texts: ArcGIS Query Layer performance documentation or feedback

texts:
I first tried asking this question in the ArcGIS forums but as is common I got a total of 0 responses. I'm willing to accept this may be due to a lack of knowledge on the subject, but I'm optimistic. I will cross-pollinate responses or just delete the ArcGIS post if it remains untouched.

I've just started using Query Layers for read-only maps served through ArcGIS Server (10 SP2) reading data from Oracle 10g. I love it! I've always felt oddly uneasy about SDE and this seems a good sidestepping opportunity.

In my tests performance looks good but I would like to see some documentation (ideally) or other user feedback on Query Layer performance on live services or large datasets and how it compares to other DBMS options (SDE, SDE DC).

The only docs I've seen make a big deal out of the data being queried each time a user action occurs, but I would expect this to happen in all types of connections. Does SDE do some magical caching that a Query Layer doesn't? Does it better handle concurrent users in some way?

If anything I would expect Query Layers to improve performance at runtime (lower overhead) and reduce development / setup time (fewer SDE configuration headaches).

Any experience or informed thoughts on this would be very much appreciated. Please assume I am not concerned with versioning / multi-user editing, or any of SDE's other non performance-related benefits.


texts: UTFGrid and line features: UTFGrid not consistently triggering mouse events for line features

texts:
[Update: My first solution works (notes below), but I'm still interested in responses to the second solution I suggested.]

Original Post:

Has anyone experienced and resolved difficulties with UTFGrid tiles of linestring feature data?

I'm encountering three issues:

The "mouse targets" for the line features are extraordinarily thin, and I have to circle the mouse around the lines to trigger anything at all. Even when I find the target area, it's difficult to follow the line and stay on target because the hit area is razor-thin.
In some places, the linear grid target seems to be a slightly offset from the rendered line feature. The trouble here is it makes the interaction counter-intuitive for end users.
Line features seem to be entirely missing from some areas of the grid. If I slowly move the mouse pointer in a circle over a line, I can usually catch the hit area and get a hover popup. However in other areas--where the features are not even dense--there is simply no hit response, leaving me to believe the feature is entirely absent from the grid.

I'm using Tilestache and Mapnik to render my grid tiles, and these are the corresponding mapnik stylesheet (xml) and tilestache config (cfg/json) file entries, respectively:

Mapnik stylesheet for linear UTFGrid layer (XML):

    <Layer name="watermains_utfgrid" srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs +over">
        <StyleName>style</StyleName>
        <Datasource>
          <Parameter name="type">postgis</Parameter>
          <Parameter name="host">localhost</Parameter>
          <Parameter name="port">5432</Parameter>
          <Parameter name="dbname">db_name</Parameter>
          <Parameter name="user">db_user</Parameter>
          <Parameter name="password">db_password</Parameter>
          <Parameter name="table">watermains</Parameter>
        </Datasource>
    </Layer>

Tilestache config entry for UTFGrid layer (JSON):

"layers":
  {
    "watermains-utfgrid":
    {
        "provider":
        {
            "name": "mapnik grid",
            "mapfile": "/etc/TileStache/map_styles/utfgrid_layers.xml",
            "fields": ["handle","watermaindia","type","est_year","subdivision"],
            "layer_index": 0,
            "scale": 4
        }
    }
 }

I have two ideas for how I might resolve these issues:

Create a view of the layer table that buffers the geometry. Rendering it as a polygon should yield a larger hit area, then I can reference the view for the Mapnik grid layer.

[Update per above] Creating a view that buffers the geometry (I used a buffer radius of "4", which is 4 meters since the data is in EPSG3857) solved all three problems, however the grids rendered considerably slow and I suspect at an undesirable expense of CPU/RAM. Storing the view as a table made a remarkable improvement to rendering speed and resource consumption. Although IMO it seems unnecessary to do this for every dataset of line features. As such I'm still hoping someone may be aware of something I should be doing differently, such as adjusting the scale value as proposed in my second idea, below.

"Tweak" the scale value in the tilestache cfg entry. This one is more mysterious, as I don't understand how that scale value affects the grid. For this reason, I actually suspect tweaking the scale is the way to go, and it's on this point I'm optimistic someone might have some useful insight.

I'm pretty sure option #1 will provide a work-around, but before I setup a bunch of grid layers against views with buffered geometries, I wanted to see if anyone else has bumped into similar issues and resolved it by adjusting the tilestache scale value ..or anything else I may have overlooked!


texts: Forcing Hash Join Query Plans using PostGIS

texts:
I'm currently testing PostGIS 2.0 in a Greenplum Appliance (Postgres 8.2 based) and have hit an incompatibility between PostGIS and the MPP Greenplum environment.

The basic query below can't be processed in parallel because it's query plan uses a nested loop for the ST_Contains operator. Using a hash join would solve the problem.

SELECT pnts.gid, bdys.pid
  FROM bdys
  INNER JOIN pnts
  ON ST_Contains(bdys.geom, pnts.geom);

I've tried using the following SET commands unsuccessfully

SET enable_nestloop = OFF;
SET enable_hashjoin = ON;

Is there any way to force the query or a variant of a point in polygon query to use a hash join?


texts: Finding pixel location in raster using coordinates?

texts:
I'm not a GIS person. I haven't found or at least identified a solution to my problem, so decided to ask.

I have a big rasters and a list of thousands of coordinates. I need to find the pixel location corresponding to a given lat/lon pair in the raster, and then clip an NxN rectangle with the pixel located in the middle of it and save the clip on disk. I haven't for the life of me figured out yet how I could find the pixel location using coordinates, and I've been on it a few days now. I found that with snappy I should be able to do something like that, but unfortunately it doesn't run on my system (anaconda3 with Python 3.6).

How can I go about solving this problem?


texts: Clipping vs Filtering Images with a Polygon - Google Earth Engine

texts:
I need to analyze the area defined by a polygon that I drew (call it 'X') of the images of Modis Terrra. What is the difference between filtering bounds and clipping to define the area of interest?

I'm confused about the following alternatives:

option 1.

var areaofinterest=ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI');
.filterBounds(X) '

or
option 2.

var aux=ee.ImageCollection('NASA/OCEANDATA/MODIS-Aqua/L3SMI')';
var areaofinterest=aux.clip(X)  ;

The main question is: what are the differences between the two methods?

When I do option 1 and then Map.AddLayer(areaofinterest) I don't get displayed the area of interest, but the map of the product over the whole world - so I guess something is wrong.
When I do option 2 I get the error: aux.clip is not a function
in <global>


texts: Finding Height and Width of Polygons using ArcGIS Desktop?

texts:
I am using ArcGIS 10.4.1

I have been asked by my employer to take a shapefile and produce a table that shows the height and width (in meters) of each polygon in a shapefile. I have explored the calculate geometry option and the only options it produces is area, perimeter and coordinates of centroid. I have also tried using the Minimum Bounding Geometry tool but this does not produce the result I want either.

Below is an example image of the shapefile. To be clear, I am interested in finding the height and width of each individual polygon.


texts: Plugin widget, include file extension by default

texts:
I've got a Plugin widget box I'm trying to configure. In one step the user defines the path to where they will export the final layout as a .pdf file. I'm using the following code to configure the box:

class CreatePrintLayoutAndExportMap(QgsProcessingAlgorithm):
    OUTPUT_PDF_PATH = 'OUTPUT_PDF_PATH'
    OUTPUT_PDF_FILENAME_DEFAULT = 'Choose file path for pdf output'

    def initAlgorithm(self, config):

        self.addParameter(
            QgsProcessingParameterFileDestination(
                self.OUTPUT_PDF_PATH,
                self.tr(self.OUTPUT_PDF_FILENAME_DEFAULT),
                ".pdf"))

Here is what the box looks like as a GUI tool:

This is good. The issue here is when a user specifies the destination, the file extension is not included and they must type "...filename.pdf" and the file extension is not manually included. Here is what it looks like in action:

and then back in the widget box. Notice the .pdf file extension is not included automatically:

If I run the tool, the file is saved as some weird file with no file type. If I explicitly include .pdf at the end of the filename (ex: "...filename.pdf"), I get a .pdf file

I've been looking at documentation for QgsProcessingParameterFileDestination as well as for QgsProcessingParameterRasterDestination as a comparison. The interesting thing is that if I export a raster file using the RasterDestination class, the .tif file extension is automatically included, removing confusion for the user.

Here is an example of a Raster output. I did not have to manually type '.tif' for it to be included:

Is there a way to automatically include a file extension when using QgsProcessingParameterFileDestination class?


texts: Unable to use ee.Kernel Google Earth Engine pyhon API

texts:
I'm trying to use the ee.Kernel Google Earth Engine pyhon API in the following code:

import ee
import ee.mapclient

ee.Initialize()
ee.mapclient.centerMap(15.88, 41.78, 10)

def maskS2clouds(collection):
    qa = collection.select('QA60')
   collection.updateMask(sombra.eq(0)).updateMask(nubes.eq(0).updateMask(cirrus_detected.eq(0)))
    cloudBitMask = 1 << 10
    cirrusBitMask = 1 << 11
    mask = qa.bitwiseAnd(cloudBitMask).eq(0)and(qa.bitwiseAnd(cirrusBitMask).eq(0))
    return collection.updateMask(mask).divide(10000)

polygon = ee.Geometry.Polygon([[
[16.440171600655958,42.033470726105556],
[15.511826874093458,42.12316884368949],
[14.621934295968458,41.689821328579825],
[15.441789032296583,41.00320702735848],
[16.246537567452833,41.36906568596335],
[16.440171600655958,42.033470726105556]]])

collection = ee.ImageCollection('COPERNICUS/S2').filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)).filterDate('2018-01-1','2018-06-30').map(maskS2clouds).filterBounds(polygon)

medianpixel = collection.median()
medianpixelclipped = medianpixel.clip(polygon)#.divide(1000)

ee.mapclient.addToMap(medianpixelclipped, {'bands': ['B4', 'B3', 'B2'], 'min': 0.0, 'max': 0.3, 'gamma': 1.5}, 'Sentinel_2 mosaic')

ndwi2 =  collection.mean().normalizedDifference(['B3','B8'])
ndwi2Viz = {'min': -0.4, 'max': 0.0, 'palette': ['FF0000', '00FF00']}
ndwi2masked = ndwi2.updateMask(ndwi2.gt(-0.125))
pixelcountndwi2 = ndwi2masked.int().connectedPixelCount(256)
ee.mapclient.addToMap(pixelcountndwi2, {'palette': ['FF0000']}, 'Sentinel_2 mosaic')
patchidndwi2 = ndwi2masked.connectedComponents(ee.Kernel.plus(1), 256)

The problem is that this error occurs:

  Undefined variable from import: Kernel

I didn't find any solution. How can I resolve this problem?


texts: How to automate downloading "slices" of data from Web Coverage Service (WCS)?

texts:
I'd like to automate downloading "slices" of data from a WCS service. This is a WCS 1.1.2 service, supporting netCDF. The data are updated daily - forecasts from an atmospheric regional model - and I'd like to automatically download netCDF "slices" of the data for a fixed geographical bounding box and for a couple of days out into the future. I am quite happy processing the netCDF once I've downloaded it.

The obvious brute-force approach is to read the OGC standard and devise the HTTP request myself, and code it up in some scripting language. I suspect I'll end up doing that. But it just feels wrong.

Can anyone recommend a WCS client library that I can use from code? Or a WCS client that I can run from a command-line? I know I could probably script QGIS or Grass or ArcGIS or something, but that seems too heavy a requirement for this task.

Can GDAL do it? Google has turned up the BADC Python WCS library (PDF) (now in OWSlib?) which I will take a look at. Java and GeoTools might be a possibility. And there's the OGC list of WCS clients which I could trudge through.

This is on the Windows platform. A .Net library would be ideal :)

I haven't linked to the actual WCS service as it's not yet publicly operational. But here is an edited GetCapabilities request below, in case anyone's interested in the details

<Capabilities xmlns:xlink="http://www.w3.org/1999/xlink"
  xmlns:ows="http://www.opengis.net/ows/1.1" xmlns="http://www.opengis.net/wcs/1.1"
 version="1.1.2">
<Contents>
<CoverageSummary>
<ows:WGS84BoundingBox crs="urn:ogc:def:crs:OGC:2:84">
<ows:LowerCorner>+0.0000 -90.0000</ows:LowerCorner>
<ows:UpperCorner>+358.8750 +90.0000</ows:UpperCorner>
</ows:WGS84BoundingBox>
<SupportedCRS>urn:ogc:def:crs:EPSG::4326</SupportedCRS>
<SupportedCRS>urn:ogc:def:crs:OGC:2:84</SupportedCRS>
<SupportedFormat>image/netcdf</SupportedFormat>
<SupportedFormat>application/x-netcdf</SupportedFormat>
<Identifier>------------</Identifier>
</CoverageSummary>


texts: Extracting metadata from Sentinel-2 images and specific dates

texts:
I am looking for a way to know how many images I have in a specific time frame.
(I will later want to have all of their indexes and dates written in a list, so I can address them later on).

I have filtered the ImageCollection to 1 month, which should give me approx. 6 images for Sentinel-2 (revisit time should be 5 days), but I get a ridiculously high number.
Also, I am not able to extract the dates as I need (first and last date of the collection).

var ndvi_s2 = ee.ImageCollection('COPERNICUS/S2_SR')
  .filterDate('2019-01-01', '2019-02-01')
  .map(function (image) {
    return image
      .clipToCollection(region)
      .updateMask(image.select('QA60').not())
      .normalizedDifference(['B8', 'B4'])
  })

// Limit the collection to the 10 most recent images.
var recent = ndvi_s2.sort('system:time_start', false).limit(5);
print('Recent images: ', recent);

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

// Get the date range of images in the collection.
var range = recent.reduceColumns(ee.Reducer.minMax(), ["system:time_start"])
print('Date range: ', ee.Date(range.get('min')), ee.Date(range.get('max')))

// Sort by a cloud cover property, get the least cloudy image.
var image = ee.Image(recent.sort('CLOUD_COVER').first());
print('Least cloudy image: ', image, ' Cloud Cover: ', image.get('CLOUD_COVER'));

These are the results:

Any ideas? Why do the "ALL" numbers is so big? And why can't I get the dates?

By the way - is there a way of seeing the NDVI values for the images I have printed out? I mean - see it in the console.


texts: Filltering 2nd combobox according to 1st combobox using PyQGIS?

texts:
I have a layer of districts. Each district belongs to a different state. I want to write code that once you choose a state in the 2end combobox, you can only see districts from that certain state in the 3rd combobox.
I wrote the following code:

       def myFunction(self, layer, featureSet2):
    Feature_list3 = []
    del Feature_list3[0:len(Feature_list3)]
    indx = self.dlg.comboBox.currentIndex()
    Data = featureSet2[indx]
    layer.removeSelection()

    # Fill 3rd combobox
    self.dlg.comboBox_3.clear()
    for feature in layer.getFeatures():
        name = feature["NAME_2"]
        if Data==feature["TYPE_2"]:
            #print (Data)
            Feature_list3.append(name)

    self.dlg.comboBox_3.addItems(Feature_list3)
    # Fill 3rd combobox

def run(self):
    """Run method that performs all the real work"""
    layer = iface.addVectorLayer("C:/Users/carmel.han/Documents/QGIS/qgis_sample_data/shapefiles/regions.shp",
                                 "regions", "ogr")

    Feature_list = []
    Feature_list2 = []

    selection = []

#Fill 1st combobox
    self.dlg.comboBox.clear()
    for feature in layer.getFeatures():
        name = feature["name_1"]

        Feature_list.append(name)
    featureSet = set(Feature_list)
    self.dlg.comboBox.addItems(featureSet)
# Fill 1st combobox

# Fill 2end combobox
    Feature_list2 = []
    self.dlg.comboBox_2.clear()
    for feature in layer.getFeatures():
        name = feature["TYPE_2"]

        Feature_list2.append(name)
    featureSet2 = set(Feature_list2)
    self.dlg.comboBox_2.addItems(featureSet2)
    Feature_set2list2 = list(featureSet2)
# Fill 2end combobox

    slotLambda = lambda: self.myFunction(layer, Feature_set2list2)
    self.dlg.comboBox_2.currentIndexChanged.connect(slotLambda)

The if statment does not seem to work,it only gives me the districts of the first state, and does not change after further changes. can somebody point out the problem or show me how it's done?


texts: Removing same name layers from layer table in JMapframe in geotool?

texts:
I am making a geotool application in which I read the data from postgresql database. In this application if i make a buffer on a selected polygon then a buffer is created and buffer layer shown in layer table of JMapFrame. After that i make another buffer on another polygon then there are two layers with same name shown in the layer table in JMapFrame and both the layers shows the buffer on last selected polygon. How can I replace the same name layer with the new one from the layer table in JMapFrame.

My code runs on the click of "Buffer" Button

And Result on Jmapfrmae

In this JMapFrame all the khasra_buffer layers are the same i.e. the current buffer which is created at Id 80, Thus replacing previously created two buffers.

I want that old khasra_buffer layer should be removed and the only buffer layer that should be shown is of the current buffer on JMapframe.


texts: Adding more fields in locator style to creat customized address locator

texts: I have a database following the structure of "US Dual Range Addresses" and I want to, in addition of the fields used in dual range addresses, add more fields and I want to know what are the changes I have to make or the steps I should follow.


texts:
If I start a project and want to use qgis2web:
Nothing happens, I see only the hourglass. Then I must restart qgis.

If I start only qgis and then start qgis2web, after click to "leaflet" I get an error message:

qgis2web produced an error:

Traceback (most recent call last):
File "C:/Users/juko/.qgis2/python/plugins\qgis2web\maindialog.py", line 318, in previewMap
preview_file = self.createPreview()
File "C:/Users/juko/.qgis2/python/plugins\qgis2web\maindialog.py", line 280, in createPreview
dest_folder=utils.tempFolder()).index_file
File "C:/Users/juko/.qgis2/python/plugins\qgis2web\leafletWriter.py", line 98, in write
folder=dest_folder)
File "C:/Users/juko/.qgis2/python/plugins\qgis2web\leafletWriter.py", line 296, in writeLeaflet
if cluster[count]:
UnboundLocalError: local variable 'count' referenced before assignment

I tried to delete and reinstall the plugin with the qgis installer: same error.
I tried to find a solution within the forum, but it is difficult as a not-programmer for me to see, if some of the error-solutions could help me.

Any idea?
QGIS under Windows 10


texts: QGIS 3 Geopackage - different attributes to same polygon

texts: I've got a geopackage with one set of polygons and an attribute table with lots of fields  - I have calculated a lot of them with reference to each other. They concern three different types of commercial property: retail, industrial and office. As they are, the attribute table is a pain to use. I'd like to split it up into 3 different objects. I could do this just by creating 3 identical polygon datasets with 3 different attribute tables, but this seems inefficient. Is there a way that I can save the data so that a single polygon dataset can relate to 3 different attribute tables - similar to structuring related tables in a database?


texts:
I'm using QGIS 1.9.0-Master, code revision 1570240.

I know I shouldn't be using it for 'production' but I'm working on a map layout, and all of the new composer features are very helpful!  As of the version available in OSGEO4W (I think) last Thursday, I was able to export my map properly to either PDF or an image. However, when I try with a version from yesterday or today, the font is not being exported properly. See examples below of proper exporting (top image) and the current export (bottom image):

I tried via OSGEO4W to roll back to a previous version, but it did not install properly for some reason.

Thoughts?


texts: Right-Align Legend Feature Counts in ArcGIS Pro

texts:
I would like to align the feature counts in my legend, to look like this legend, generated in ArcGIS Online:

which I believe is much easier to read than this one, generated in ArcPro:

Is there any way in ArcPro to get the feature counts right-justified, like in the top image?

On a separate note, is anyone aware of a way to adjust the font size of the heading (in this case, "Structural Concerns") without also adjusting the size of all the other legend labels?


texts: Empty signature file after i.gensig in GRASS

texts:
I'm trying to make a supervised classification in Grass 6.4.3 under Linuxmint 16.

After digitizing the training areas, I transformed them from vector to raster.
Then I've done the following:
just in case:
g.region rast=train_area_2013_03_05_raster@PERMANENT,2013_03_05_toar4@PERMANENT

then

i.group group=group subgroup=sub_group input=2013_03_05_toar4@PERMANENT,2013_03_05_toar5@PERMANENT,2013_03_05_toar3@PERMANENT

and

i.gensig trainingmap=train_area_2013_03_05_raster@PERMANENT group=group@PERMANENT subgroup=sub_group signaturefile=superv_class

but the resulting signature file, only contains one character: "#".

The output of i.gensig in the console is the following:

Finding training classes...
3 classes found
Calculating class means...
Calculating class covariance matrices...
Signature 1 not invertible
Signature 2 not invertible
Signature 3 not invertible
i.gensig complete.

Obviously if I try:
i.maxlik group=group@PERMANENT subgroup=sub_group sigfile=superv_class class=results_class

I obtain the following error message:

ERROR: Unable to read signature file <superv_class>

I've been searching, but I couldn't find the solution.

I've tried: this Problem running i.gensigset in Grass. Any ideas?

My problem is similar to this: i.maxlik cannot read i.class output signature file. But here, there's no solution suggested. As in this case, I also used i.cluster to generate an unsupervised signature file and i.maxlik reads this perfectly well.

Here there's another similar question, but again, with no solution...

What does Signature X not invertible means?

Any idea?
Thanks for your help!


texts: Calculate shortest path and distance between multiple OD pairs using pgr_dijkstra one-to-one

texts:
My goal is to find the distance and shortest path between 1000+ Origin-Destination pairs in a road network. I have a table with the pair ID, the origin vertex ID and the destination vertex ID for each pair.  Everything works fine when I use the pgr_dijkstra, where I get a path and a agg_cost for the following example:

SELECT * FROM pgr_dijkstra(
'SELECT gid as id, source, target, length as cost FROM route1',
48302, 44382, FALSE);

What I would like to do is replace the origin and destination vertices (48302, 44382) by the Origin and Destination vertex id columns from my table, matched by the pair ID.  In other words, I'd like to get the distance and path for all my OD pairs in a single query.


texts: How to call a temporary layer from another python script in QGIS?

texts:
I wrote several python scripts called by my main.py script.
In one of those scripts I create a vector layer by polygonizing a raster layer. As the layer is really heavy and as I don't need it as an output, I create it in memory.

result = processing.runalg('gdalogr:polygonize',layer,"lu_majorit", None)
polygonLayer = processing.getObject(result['OUTPUT'])

Then I can use the layer in my file. But I would like to use this layer from another script and I don't succeed.

Here is what I've tried :

in my first script :

result = processing.runalg('gdalogr:polygonize',layer,"lu_majorit", None)
polygonLayer = processing.getObject(result['OUTPUT'])
return polygonLayer

in my main:

polygonLayer=script1(...)
script2(polygonLayer,...)

in my second script :

cLayer =QgsVectorLayer("Polygon?crs=epsg:4326", "OUTPUT", "memory")

or:

cLayer =QgsVectorLayer("polygonLayer", "OUTPUT", "ogr")

Those two ways failed to load the layer.

I tried also by giving the path with :
script 1 :

myfilepath= polygonLayer.dataProvider().dataSourceUri()
return myfilepath

script2:

cLayer =QgsVectorLayer(myfilepath, "OUTPUT", "ogr")

This way gives me an invalid syntax.

Does someone know how to do ?


texts: Qgis 3.4 - How could I add a QDockWidget to Qgis interface using Python

texts:
I was looking for a way to add a panel and I found QDockWidget but I don't seem to find the way to add a new panel since examples I've seen use a Qgis tool to create these panels or use QDockWidgetto show or hide existent panels by their name.

So the question is: How could I add a custom QDockWidget as a Qgis Panel with the content I tell it to load?

This is the way I've been able to create it, it can only be appended to the plugin's interface.

I want it to be placed in the Qgis interface like this one:


texts: Using PostgreSQL trigger on table that updates field?

texts:
Depending on values on 2 fields A and B, I would like a field C to be updated or filled with a specific value.

For example :

create or replace function periode() returns trigger as
$BODY$
begin
IF A >= x AND B <= y
THEN insert into C values ('z')
END IF;
RETURN NULL ;
END ;
$BODY$

LANGUAGE plpgsql;
DROP TRIGGER IF EXISTS trg_periode ON temp ;
CREATE TRIGGER trg_periode AFTER INSERT or UPDATE on temp
FOR EACH ROW EXECUTE PROCEDURE periode();

It does not work : 'ERROR : column A doesn't exist'. I can't figure out why this message appears.

Maybe there is a better way to create this trigger?


texts:
I was given an RGB GeoTIFF orthomosaic image that is just over 10GB. I pulled it into QGIS and tried clipping it to get under 10GB (to use in another program that only allows up to 10GB).
Though when I clip it or do any other process, the output file size is over double of the orginal (24GB). I tried setting the zero values to null, but that has minimal effect. I need to keep the image as an RGB, so I cannot export it as a jpeg or png.
Has anyone else run into this problem and solved it, if yes then how?


texts: Obtaining closest distance using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.5.

I have a shapefile containing say 10 locations from a country with latitude and longitude.

How can I look at each location (point) and getting the distance from each of the other 9 locations (points) and storing the "closest distance" and the location that appears the closest?

In the end, I would have a table showing the 10 locations along with the closest location name with the distance.


texts:
I have a OSM server running (Geoserver + PostGIS + OSM Planet Data) on my own server.

Now, I need to offer to my users the same functionality of openstreetmap.org to extract data from the database using a web interface.

How can I do this?


texts: Unable to connect to ArcGIS Server in ArcGIS Pro

texts:
I have an ArcGIS Server URL that I connect to no problem with ArcMap, but when I try to add it in ArcGIS Pro using any of the connection means available--it fails to connect. It asks me for my credentials but gives me an "unable to connect" error.
Has anyone else had this issue when migrating from ArcMap to ArcGIS Pro?


texts: OpenLayers + MapServer - Finding reason for random missing tile

texts:
I am testing an application this is developed in OpenLayers 2.13 with Chrome, and I am using MapServer as the map provider.

I have noticed some tiles occasionally does not render, when I pan/zoom the map. There is no specific way to recreate this problem on my application, the missing tiles are random. But I can say that the problem is common enough for most users to notice.

MapServer is running on an dedicated core i7 desktop with 16GB of RAM. I do not understand if it is a hardware issue, software issue, or load handling issue.

I have found some links that specifies the solution to hide the pink image
XYZ Layer in OpenLayers has missing tiles

But I am interested to know what is the cause of missing tiles? What are the possible factors causing missing tiles.

UPDATE

Sometimes it will also appear that the tiles are loaded or what appear to be loaded correctly, turned from map tiles into blank tiles.