texts: Keeping output from ModelBuilder in projected coordinate system of input?

texts:
My model input is a point feature class with groups of points which need to be sorted according to their unit_id and then a Minimum_Bounding_Geometry tool run to generate a polygon for each unit_id which is then saved in a Feature Dataset.

This all works in the model attached except my input point layer is a projected UTM coordinate system. However the output polygons are not projected. The reason this matters is that I want to use the auto-generated shape_Area attribute value but for that value to be a sqm area, I believe I need it to be in a projected coordinate system.


texts: arcpy connect to accdb with 64bit background processing

texts:
I'm trying to read from a accdb MS Access Database via python in a ArcGIS Toolbox tool. I already established a database connection (via Add OLE DB in ArcCatalog) which I can use in a 32bit python process. I have the 64bit background processing installed and don't want to set

self.canRunInBackground = True

in the Python Toolbox.
When writing something like

with arcpy.da.SearchCursor(test.odc.table_name,[field_name]) as sc_rows:
    for sc_row in sc_rows:
        print sc_row[0]

it works perfectly fine in a 32bit python environment but in 64bit python I get the error message:
RuntimeError: cannot open

I have Windows7 64bit, ArcGIS 10.4.1 including the 64bit background processing but only a 32bit MS Office. Might that be the problem?

Is there any  way to access MS Access accdb databases from python having the 64bit background processing installed?

The documentation of ArcGIS 10.1 states that ole-db connections are not suported in 64bit background processing (http://resources.arcgis.com/en/help/main/10.1/index.html#/Background_Geoprocessing_64_bit/002100000040000000/) whereas in the documentation of ArcGIS 10.4 it is not listed as unsupported data type (http://desktop.arcgis.com/en/arcmap/10.4/analyze/executing-tools/64bit-background.htm).


texts: Topographic Wetness Index difference between SAGA and TauDEM

texts:
I calculated Topographic Wetness Index (TWI) using SAGA and TauDEM. Overall trend is similar but SAGA returned larger values compared with TauDEM. Sometimes the deviations are as large as 7 or 8. I am using 1 m resolution DEM. I suppose the calculation methods are the same.

What is causing these differences?


texts: Leaflet.draw: add attributes and save to file

texts:
I want to create a simple web page, where user can create features (with Leaflet.Draw), add attributes (ex. name/description) and save created features to file, ex. geojson. So far creating, adding attributes and saving is working, but saved geojson files does not have added attributes. Maybe they are not added or mayby I need to update drawnItems in some way...
I have learned from these resources:

Adding comments using leaflet draw

Leaflet Draw - add title to marker or shape?

Here is my code:

Export Features

//Leaflet.draw
var drawnItems = L.featureGroup().addTo(map);
var drawControl = new L.Control.Draw({
 draw: {
  circle: false,
  rectangle: false,
  polyline: false,
  polygon: false
 },
 edit: {featureGroup: drawnItems}
}).addTo(map);

map.on('draw:created', function(e) {
  var type = e.layerType;
  var layer = e.layer;
  var idIW = L.popup();
  var content = 'Shape Name
<input id="shapeName" type="text"/>

Shape Description
<textarea id="shapeDesc" cols="25" rows="5"></textarea>

<input type="button" id="okBtn" value="Save" onclick="saveIdIW()"/>';
  idIW.setContent(content);
  idIW.setLatLng(layer.getLatLng());
  idIW.openOn(map);
  drawnItems.addLayer(layer)
});

function saveIdIW() {
  var sName = $('#shapeName').val();
  var sDesc = $('#shapeDesc').val();
  var drawings = drawnItems.getLayers(); //drawnItems is a container for the drawn objects
  drawings[drawings.length - 1].title = sName;
  drawings[drawings.length - 1].content = sDesc;
  map.closePopup();
};

//Export
document.getElementById('export').onclick = function(e) {
  // Extract GeoJson from featureGroup
  var data = drawnItems.toGeoJSON();
  // Stringify the GeoJson
  var convertedData = 'text/json;charset=utf-8,' + encodeURIComponent(JSON.stringify(data));
  // Create export
  document.getElementById('export').setAttribute('href', 'data:' + convertedData);
  document.getElementById('export').setAttribute('download', 'drawnItems.geojson');
}

I have created a http://jsfiddle.net/falqn/svha4bpt/ but I do not know why 'save' button do nothing and there is an error: 'Uncaught ReferenceError: saveIdIW is not defined'. I'm new to jsfiddle, maybe I'm doing something wrong. On my test page (local) 'save' button is closing the popup and everything seems to be OK, but in saved geojson file there is no 'Name' and 'Description' attributes in it.
I think I'm very close to really easy way to create/edit/save features with Leaflet.Draw. Please help.


texts: Can QGIS for Windows use ticks instead of crosses in the Table of Contents?

texts:
In a fairly old blog, Martin Dobias shared a post detailing the use of the QGIS Layer Tree API where one of the images shows the Table of Contents (or Layers Panel in QGIS jargon) using ticks to show if a layer is enabled or not:

This seems to be available in QGIS for Linux users but for Windows, it uses crosses:

Is it possible to use ticks for a Windows version of QGIS?


texts: Create a line chart from a dictionary of values in google earth engine

texts:
I am currently trying to create a line chart for each land-cover type across three bands I selected from MODIS.

Problem is most of the ui.Chart functions require a feature or an image. Please help?

// Create clipping function
var clipper = function(image){
  return image.clip(geometry);
};

// Load the MODIS data that is being processed

// Clip the image collection
var filtered = ee.ImageCollection('MODIS/006/MOD09A1').map(clipper)
                  // Sort dates in descending order
                  .sort('system:time_start',false);
print(filtered);

// Create a visualization color scheme
var medVis = {
  min: 0,
  max: 8000,
  // You have to specify one band to use a palette
  bands:['sur_refl_b04'],
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

// Load the first image
Map.centerObject(geometry,8);
Map.addLayer(filtered.first(),medVis,"Loading the most recent Image");
Map.addLayer(geometry, {}, 'Shapefile Outline');

// Get the recent image alone
var singleImage = filtered.first();
print(singleImage);

// Merge the featureCollections that have been created
// Create your own polygons here to make the code reproducible.
// It can be polygons or points that refer to a particular land cover type
var training = forest.merge(water).merge(bare);
print(training);

var mean = singleImage.select(['sur_refl_b02','sur_refl_b03','sur_refl_b04'])
                      .reduceRegions({
                        collection: training,
                        reducer: ee.Reducer.mean(),
                        scale: 500
                      })
                      .reduceColumns({
                        selectors: ['sur_refl_b02','sur_refl_b03','sur_refl_b04','class'],
                        reducer: ee.Reducer.mean().repeat(3).group({
                              groupField: 3,
                              groupName: 'class'})
                        });

print(mean, 'Band Means');

Is there also any way in which I can retain the names of the selectors in the dictionary created as well?

I have added a link to code based on Rodrigo's comment: https://code.earthengine.google.com/a84485d5cba393122e2a1c4c5ae17690


texts: Creating a polygon that shows the specific elevation above a river using a DEM

texts:
My goal is to determine the possible area covered by interpreted Tertiary drainages.  To keep my example simple, I have a DEM and a River polyline.  I would like to see the area a polygon (or raster) would cover if it was 100 feet above my River polyline.  Being that the river decreases in elevation as it drains, the polygon would cover most of the valley but some might spill over the top and could possibly show where old drainages might have occurred.  I have attempted to google this but have not found anything.  Seems simple enough to do.

Can someone lead me in the right direction?

I am using ArcMap 10.5


texts: Problem with Attribute Assistant Expression Command

texts:
I am using the Attribute Assistant add-in to perform multiple steps when a new feature is created. Two of the important steps are: generate a unique, incremental ID (project ID) and then use that project ID to build a larger docket number.

I have gotten the generate ID field to work properly and create a number in a field called PROJECTID. I then also have a field called COMPLAINT which tracks the type of complaint. The expression I want to build is [COMPLAINT] & "-2016-" & [PROJECTID], so the docket number will look something like VTW-2016-1.

I can get the the expression to work when I just use [COMPLAINT] & "-2016-", but when I add on the [PROJECTID] to the end of the expression, the tool fails and the field remains Null.

I will attach two images below, one of my DynamicValue table and one of the log that is tracking the error.

We are using ArcGIS 10.2.2 in an SDE database on SQL server.


texts: postgresql/postgis division by zero error

texts:
I have a query where I am calculating acreage of two intersecting tables and what percentage that new acreage is compared to the old acreage.here is the error

ERROR:  division by zero
********** Error **********

ERROR: division by zero
SQL state: 22012

query

select p.objectid,acreage old_acreage,p.shape,sum(st_area(st_intersection(p.shape,lulc.shape))/43560) as new_acreage,
    sum(st_area(st_intersection(p.shape,lulc.shape))/43560)/acreage as agr_use_percentage

from base_layers.modiv_parcels_2015 p join base_layers.lulc_huc_hl_clip lulc on st_intersects(p.shape,lulc.shape)

where lu12 in(2100,2200,2300,2400,2140,2150) group by p.objectid,acreage,p.shape

what can be done to circumvent this error?


texts: How to do a looping to Interpolate a huge Rainfall Time-series in GRASS GIS or R

texts:
In GRASS GIS I have access to a table in PostgreSQL Database with this structure:

Table Name: rainfall_time_series

Columns: id | year | month | day | value_mm | x | y | gauge_code

Setting g.region to the Watershed limits, I need create for each day one interpolated raster, representing the rainfall. This times-series is very huge, has 17 years, my question is: is it possible to do this operation with a looping? it is probably easier in R, but I'm newbie in this software.

Thanks in advance


texts: Calculating right area in polygons using DEM and ArcGIS Desktop?

texts:
I have used infra red photos to identify territories in Norway where I have drawn polygons for different vegetation. I wonder if there is a way to use the polygons with elevation data to get the exact area since the territories varies a lot from high to low.

I also wonder if there is a easier way to identify the vegetation and to get the exact area with help of elevation data.

I´m using Windows 8.1 and I´m working in ArcGIS.


texts:
I am trying to show a list as a datagrid full of strings in a dockable window.

I can bind the list to the datagrid but I can't for the life of me get the datagrid to show in the dockable window.

I have the datagrid modifier as public

Any ideas?

public class NotLabelled : ESRI.ArcGIS.Desktop.AddIns.Button
{
    public static IDockableWindow s_dockWindow;

    public NotLabelled()
    {
    }

    protected override void OnClick()
    {
        UID dockWinId = new UIDClass();

        dockWinId.Value = ThisAddIn.IDs.ResultsWindow;
        s_dockWindow = ArcMap.DockableWindowManager.GetDockableWindow(dockWinId);

        s_dockWindow.Show(!s_dockWindow.IsVisible());

        ResultsWindow.CheckForRoadsNotLabelled();
    }

    protected override void OnUpdate()
    {
        Enabled = ArcMap.Application != null;
    }
}

/// <summary>
/// Designer class of the dockable window add-in. It contains user interfaces that
/// make up the dockable window.
/// </summary>
public partial class ResultsWindow : UserControl
{
    public static List<string> AnnoList = new List<string>();
    public static List<string> StreetList = new List<string>();
    public static List<string> NotLabeled = new List<string>();

    public ResultsWindow(object hook)
    {
        InitializeComponent();
        this.Hook = hook;

        // load other data
        // ie xml data
    }

    /// <summary>
    /// Host object of the dockable window
    /// </summary>
    public object Hook
    {
        get;
        set;
    }

    /// <summary>
    /// Implementation class of the dockable window add-in. It is responsible for
    /// creating and disposing the user interface class of the dockable window.
    /// </summary>
    public class AddinImpl : ESRI.ArcGIS.Desktop.AddIns.DockableWindow
    {
        public static ResultsWindow m_windowUI;

        public static List<string> NotLabeled = new List<string>();

        public AddinImpl()
        {
        }

        protected override IntPtr OnCreateChild()
        {
            m_windowUI = new ResultsWindow(this.Hook);
            return m_windowUI.Handle;
        }

        protected override void Dispose(bool disposing)
        {
            if (m_windowUI != null)
                m_windowUI.Dispose(disposing);

            base.Dispose(disposing);
        }
    }
}

public static void LoadGrid(IList<string> notlabeled)
{
    DataGridView dgResults = new DataGridView();

    dgResults.Rows.Clear();
    dgResults.Columns.Clear();
    dgResults.Refresh();

    dgResults.SelectionMode = DataGridViewSelectionMode.FullRowSelect;

    dgResults.DataSource = NotLabeled; // notlabeled is binded to dgResults

    dgResults.Show(); // doesnt work
}


texts: Work around to create feature-linked annotation from standard annotation?

texts:
I have one standard annotation feature class that I have used to label the contract names of 6 feature classes. I'm in a tough position where I've been asked to remove a few hundred records from about 2000 records spanning over the 6 different feature classes. This would create 200 orphaned annotations in my standard annotation feature class and I need to have the flexibility to bring these annotation back just in case. Since the labels are meant for a data driven map book, the labels and leader lines have been painstakingly positioned manually. I don't want to have to recreate all these labels as "feature-linked" because I would have to re-position them all.

Since the text field from the annotation layer has a corresponding text field in the feature classes, I should theoretically be able to create a link between the two with a relationship table. I have no idea how to do this or if it is even possible to retro-actively create a feature-linked annotation. Has anyone done this or have any answers on how to do this?


texts: Rasterise a vector layer but relevant column is not numeric in QGIS

texts:
The software I'm using is QGis 3.16.0.
In order for a study project, I need to rasterise a vector layer (polygons) in order to classify certain classes of a certain attribute column with the process "Reclassify by table".
The attribute COVERTYPE consists of eight classes: AGR, BAR, CON, HDW, MIX, SHB, URB and WAT.
I tried to rasterise my vector field using the SAGA Tool "Rasterise". As "Attribute" I chose COVERTYPE, which is the correct one. For my other options please look at this screenshot:

Might there be any wrong settings in my configuration? I'm not that familiar with rasterising just yet...
For my rasterised layer, I just want to have eight different classes as well which I can classify then easily with the raster calculator. How do I know now, after I successfully rasterised my vector layer, which value represents which class?
By looking at the histogramm of my rasterised layer, there should be just eight different pixel values, instead I see many different pixel values, although my attribute covertype just has eight different classes. See also:

For your information, during the SAGA Rasterise process, there is the following warning: selected attribute is not numeric; generating unique identifiers instead. See also:

Could that be the problem? Firstly, do I have to give my attributes AGR, BAR, and so on numbers? That wouldn't quite be possible, because I have around 1000 objects...
Is there any other process which is better for my certain purpose? Unfortunately, the process "Rasterise (vector to raster)" doesn't even work, there is no output map.


texts:
How can I find the inner bounds of a polygon in Leaflet? I a want to plot points randomly in polygon but when I say
Layer.getBounds()

I get the bounds and I can still plot the markers approximately with in the saure i draw except most of them fall far out the intended polygons. This is because, getBounds() gives me out bounds.
Does any one have an idea how I can get inner bounds of a polygon?


texts: PostGIS Intersection and Summarise Attributes

texts:
I am usually an ArcGIS desktop user, but keen to start using PostGIS more and I have a really big bit of processing to do. Not sure what functions to use, hopefully someone can help.

I have a polygon dataset (several million features), based on a type of landuse/ landcover classfication (20 categories). I have a number of regions in another dataset.

For each of the regions, I would like to know the area of each landcover classfication.

In ArcGIS (if it was a smaller dataset) I would imagine first adding the region to each of the polygons in the attribute table using a join. Then using "summarize" on the table by region and by landcover classification.

Not sure where to start doing this in PostGIS / SQL.

Update:

Wow thanks that has been a huge help.

It has been running a long time (44 hours!) and now I get:

NOTICE:  TopologyException: found non-noded intersection between LINESTRING (coords
edited) and LINESTRING (coords edited) at (position edited)
ERROR:  GEOS Intersection() threw an error!

********** Error **********

ERROR: GEOS Intersection() threw an error!
SQL state: XX000

I assume this is a problem in the original data - just a case of reviewing the original data or can I first check the topology some how for the whole data?
Is there something about accepting certain errors / processing tolerances?


texts: Cannot acquire lock from remote geodatabases on NAS

texts:
This is all due to the remote server being a NAS, which led me to these fixes:

ArcGIS File Geodatabase import failure on TS-453Pro
https://geonet.esri.com/message/71865#71865

which were inconsistantly successful for me (and unsafe), so I'm moving production off this server.

I believe I have a permissions issue with ArcMap (10.2.2) geodatabases (file) on a remote file server. The error message is "Cannot acquire a lock." And it crops up in a number of (failed) processes:

importing features from another gdb to the remote gdb
running dissolve on features in the remote gdb
Auto Complete Polygon/Free Hand (in Create Features)

The difficulty in tracking this down is that there is no warning there are problems with locks until near the end of the process. Or in the case of Creating Features, no other operations besides Auto Complete Polygons will fail, making Polygons, Rectangles, Ellipse, and Freehand all work fine.

The only workaround I have found is to move the gdb to my desktop or our old (XP) fileserver (which gives full permissions). Following this clue along with the error message led me to discover that the remote gdb never opens lock files for the gdb or layer files when in use, only locks for files in use by editor (ed.lock). This image is from a local copy, and shows the three types of locks:
https://drive.google.com/file/d/0B-pqTCgOoeO9M1FtcVZ6R2ltSDg/view?usp=sharing

Things which did not resolve the problem:

compact gdb
Run as Administrator
set scratch to local machine

Has anyone else had this problem and found a solution?


texts: PostGIS raster functions aspect and slope slow

texts:
I am trying to compute the raster and slope values from a DGM table (8184 x 5159) imported from tif file tiled in 500x500 tiles.

In my query in which I merge the DGM table with a vector (point type) table consisting of 3000 (or 15000) points this query takes extremely long.

Explain analyze tells me for a single point to compute aspect and slope it takes more than 3s. I set on both tables an index (gist (st_convexhull(rast)) on the DGM and a gist(geom) on the point table). Any idea how to speed up this query?

Markus

EXPLAIN ANALYZE
SELECT
  ST_Value(ST_Aspect(dem.rast, 1, '32BF'), sat.geom,true) AS aspect,
  ST_Value(ST_Slope(dem.rast, 1,'32BF'), sat.geom,true) AS slope
FROM
   safer.st_dgm_tiled dem, safer.aoi_bolz_rsat_a sat
WHERE ST_Intersects(dem.rast, sat.geom) LIMIT 1;

And this is the query plan:

  Limit  (cost=0.00..1101.35 rows=1 width=64) (actual time=3332.439..33051.935 rows=10 loops=1)
       ->  Nested Loop  (cost=0.00..1101.35 rows=1 width=64) (actual time=3332.434..33051.913 rows=10 loops=1)
           Join Filter: _st_intersects(dem.rast, sat.geom, NULL::integer)
           ->  Seq Scan on st_dgm_tiled dem  (cost=0.00..3.87 rows=187 width=32) (actual time=0.103..0.119 rows=59 loops=1)
           ->  Index Scan using aoi_bolz_rsat_a_gist on aoi_bolz_rsat_a sat  (cost=0.00..3.35 rows=1 width=32) (actual time=0.006..0.007 rows=0 loops=59)
                 Index Cond: ((dem.rast)::geometry && geom)
   Total runtime: 33065.589 ms
  (7 rows)


texts: Opening Cities Revealed data from 1990s?

texts:
I have been given a lot of aerial imagery from the mid-1990s that is in '.crv' format. The software/company at the time was called 'Cities Revealed'.

Does anybody know how I'd go about opening it?

I have: QGIS, ArcGIS, AutoCAD and FME at my disposal.


texts:
I had got this error when I try run Qwater model.

  "2019-09-15T19:47:05     WARNING    Traceback (most recent call last):
  File
  "C:/Users/OS10/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QWater\QWaterPlugin.py",
  line 688, in runEpanet
                writer = GHydraulicsInpWriter(template, self.iface)
                File "C:/Users/OS10/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QWater\GHydraulicsInpWriter.py",
  line 439, in init
                self.getVirtualLines(section)
                File "C:/Users/OS10/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\QWater\GHydraulicsInpWriter.py",
  line 321, in getVirtualLines
                length = self.getSecondVirtualLineJunction(id, elevation)

  UnboundLocalError: local variable 'elevation' referenced before
  assignment"


texts: Shapefile has feature that I cannot edit on QGIS 3.10.13 on Windows 10

texts:
I am editing a map with multiple layers. The project includes a shapefile for buildings that has numerous rectangles in black with 100% transparency as a single symbol, The layer file is at the top of the list but one other layers in the project has various polygons in  Catorgarized colors at 50% transparency. These do not let the black show through.
There are other layers with polygons of categorized colors that overlay the same black shapefile and these do not cover the black. I have tried redrawing parts of both layers but the issue persists.
I have limited knowledge of QGIS. I am trying to produce a Layout for printing . As part of this I am moving features ie Buildings on a layer and adding text to some items. I am using the vertex tool. I have done several layouts with other layers without a problem but this one layer covers the features on my building layer.

I think me saying edit is not helpful. If you look at the screen shot you can see the rectangles are showing dark green. They should be black. The rectangles are in the Building layer which is at the top of the layer panel. Why does the green feature which is in another layer cover the building layer. How can I change it?
I understand the operations available but unfortunately it does not change the map. I have deleted items on the layer and redrawn them but they are still covered by the other layer. I have shut down and reloaded twice.
You will not believe this I have just gone to the layers to do a screen dump and the buildings are now showing. I didn't change anything.


texts: SQL trigger that dissolves feature class based on multilabel fields

texts:
I am about write an SQL trigger which dissolves a feature class based on 5 fields, I used to do this using dissolve tool in the Arcmap as the following screenshots:

Now, I need this process to be automated, so when the user inserts or updates the SDE database the trigger should be able to dissolve the feature class based on the above criteria.

How can I start with this job, I used to write simple SQL triggers however that isn't sufficient to fulfill this process.


texts: How to create an ATLAS using second feature iteration on QGIS Composer

texts:
I'm looking forward to get deeper with QGIS ATLAS Composer features like we do with some reporting sofwares.

Today I try to map and iterate over parcels and display buildings attributes on the side, all buildings on a slide should be on the iterated parcel.

The idea behind could be more a generic approach to list joint subfeatures from an iterated one used for the Atlas.

For my Use case:
Buildings and parcels have a geographical link and also keys to link each other.

Update : I add a better example to be more explicit. In this configuration, we could add buildings linked pictures (from the database urls using html), and so on...


texts: How to unselect elements in PageLayoutControl (using C#)

texts:
I have a PagelayoutControl with the name: PageControl. And I added a text element into MapFrame (FocusMap) of the PagelayoutControl. If the text element and MapFrame are selected (by using the button "Select Elements" on the toolbar).

Do you know are there any ways to unselect the elements and MapFrame by using code in C# ?

Please see my screenshots (in cases: selected and unselected)

Thanks and regards,


texts: Openlayers How would i get a hover label over my polygon

texts:
make a hover/mouseover label for a polygon in Openlayers?

I made a map with draw, delete and drag polygon function.

Now I need to make a hover mouseover for this.

How to get a hover label for polygon

        var polygonLayer = new OpenLayers.Layer.Vector("Polygon Layer");

        map.addLayers([wmsLayer, polygonLayer]);

        polyOptions = {sides: 40};
        polygonControl = new OpenLayers.Control.DrawFeature(polygonLayer,
                                        OpenLayers.Handler.RegularPolygon,
                                        {handlerOptions: polyOptions} );

       var myStyles = new OpenLayers.StyleMap({
            "default": new OpenLayers.Style({
                pointRadius: "${type}", // sized according to type attribute
                fillColor: "#ffcc66",
                strokeColor: "#ff9933",
                strokeWidth: 2,
                graphicZIndex: 1
            }),
            "select": new OpenLayers.Style({
                fillColor: "#66ccff",
                strokeColor: "#3399ff",
                graphicZIndex: 2
            })
        });

         var select = new OpenLayers.Control.SelectFeature(polygonLayer, {hover: true});
        map.addControl(select);
        select.activate();


texts: Plotting and/or exporting VIIRS SDR DNB data as GeoTIFF in R?

texts:
I already asked this question @stackoverflow however, perhaps it fits better here, due to the spatial nature of VIIRS data and the hdf5-eos file format.

Im looking into how to handle .h5 (hdf-eos) VIIRS DNB raster data in R (or any other opensource software). I know that GeoTIFF products are available, for monthly composites, however I need daily data, that is why I need to work with the .h5.

My data is ordered and downloaded through https://www.class.ngdc.noaa.gov/saa/products/welcome

Preferably I want to be able to plot and create time series of the data in R, but if that is not possible, then I want to export it as GeoTIFF to handle it in ArcGIS or QGIS (which is maybe also necessary to work with it in R?).

So far I am able to read in the file with

hdf<- h5file("mypath/GDNBO-SVDNB_npp_d20171101_t0110370_e0116174_b31151_c20180208224859630066_nobc_ops.h5", mode = "r")

But I cannot figure out how to access the "radiance" band, let alone plot it. I have read the "h5" package documentation without any progress.

As another option I have tried to export the h5 file to geotiff with

sds <- get_subdatasets("subset_1_d20171101_t0110370.h5")
gdal_translate(sds[17], dst_dataset = "radiance.tif", overwrite=T) #subset 17 is the radiance band i need

however the file has lost its geolocation during the transformation, when the geotiff is opened in e.g. QGIS, it does not have any location or projection.

I understand that the geolocation is stored in separate bands (latitude, longitude) however I don't understand how I can access meta data on the corner coordinates and the projection, and then pass this on to gdal_translate() or a software like ArcGIS?

Anyone knows how to deal with these type of files in R?

I seek an example in R how to plot or export this data file.


texts: How to add GML data on arcgis javascript map application?

texts:
I am using a wfs service using ajax query. So it returns me a xml based response (GML2 format) like following. (But does not return a json.)
<wfs:FeatureCollection xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:wfs="http://www.opengis.net/wfs" xmlns:gml="http://www.opengis.net/gml" xmlns:osm="http://openstreemap.org" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" numberOfFeatures="62" timeStamp="2020-10-03T16:53:48.253Z" xsi:schemaLocation="http://www.opengis.net/wfs https://ahocevar.com/geoserver/schemas/wfs/1.1.0/wfs.xsd http://openstreemap.org https://ahocevar.com/geoserver/wfs?service=WFS&version=1.1.0&request=DescribeFeatureType&typeName=osm%3Awater_areas">
  <gml:featureMembers>
    <osm:water_areas gml:id="water_areas.23">
     <osm:the_geom>
       <gml:MultiSurface srsName="http://www.opengis.net/gml/srs/epsg.xml#3857" srsDimension="2">
        <gml:surfaceMember>
          <gml:Polygon>
           <gml:exterior>
             <gml:LinearRing>
              <gml:posList>-8894836.86 5381000.5 -8894825.07 5381032.26 -8894816.35 5381039.83 -8894774.3 5381031.02 -8894716.03 5381026.55 -8894713.34 5381014.82 -8894719.1 5381002.37 -8894742.99 5380977.7 -8894759.23 5380973.33 -8894777.12 5380973.53 -8894795.06 5380979.26 -8894836.86 5381000.5</gml:posList>
             </gml:LinearRing>
           </gml:exterior>
         </gml:Polygon>
       </gml:surfaceMember>
    </gml:MultiSurface>
   </osm:the_geom>
   <osm:osm_id>76963679</osm:osm_id>
   <osm:natural>natural</osm:natural>
   <osm:landuse>reservoir</osm:landuse>
</osm:water_areas>
<osm:water_areas gml:id="water_areas.220">
     <osm:the_geom>
        <gml:MultiSurface srsName="http://www.opengis.net/gml/srs/epsg.xml#3857" srsDimension="2">
         <gml:surfaceMember>
          <gml:Polygon>
            <gml:exterior>
               <gml:LinearRing>
                  <gml:posList>-8914042.07 5382281.26 -8914035.71 5382296.35 -8914017.03 5382316.79 -8913991.17 5382332.74 -8913989.48 5382357.25 -8913968.34 5382359.11 -8913938.22 5382340.81 -8913937.73 5382325.46 -8913987.93 5382271.05 -8913983.59 5382260.33 -8913990.68 5382253.72 -8914013.55 5382261.97 -8914042.07 5382281.26</gml:posList>
              </gml:LinearRing>
           </gml:exterior>
        </gml:Polygon>
     </gml:surfaceMember>
   </gml:MultiSurface>
</osm:the_geom>
<osm:osm_id>39458623</osm:osm_id>
<osm:natural>natural</osm:natural>
<osm:landuse>reservoir</osm:landuse>
</osm:water_areas>
....
....

So I want to add this data on map in my arcgis javascript application. How can I do it?


texts: Geoserver CQL file queryCollection limit

texts:
I'm using GeoServer 2.9.0 and I need to get some result using CQL_FILTER, intersects between polygon and line String layer, cql query is as follow:

INTERSECTS(geom, collectGeometries(queryCollection('workspace:line_string_layer', 'geom', 'IN(''line_string_layer.23'')')))

The full requested url is

http://IP:8080/MapServer/workspace/wfs?service=WFS&version=1.0.0&request=GetFeature&typeName=workspace:polygon_layer&maxFeatures=500000000&outputFormat=application/json&CQL_FILTER=INTERSECTS(geom,collectGeometries(queryCollection('workspace:line_string_layer','geom','IN(''line_string_layer.23'')')))

But the result return as follow

<?xml version="1.0" encoding="UTF-8"?>
    <ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows/1.1"
     xmlns:xs="http://www.w3.org/2001/XMLSchema"
     xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="2.0.0"
     xsi:schemaLocation="http://www.opengis.net/ows/1.1
     http://IP:8080/workspace/schemas/ows/1.1.0/owsAll.xsd">
     <ows:Exception exceptionCode="NoApplicableCode">
         <ows:ExceptionText>java.lang.IllegalStateException: The query
              in queryCollection returns too many features, the limit is
              1000 The query in queryCollection returns too many features,
              the limit is 1000
         </ows:ExceptionText>
      </ows:Exception>
</ows:ExceptionReport>

How to increase this queryCollection limit or am I doing something wrong?


texts: Difference between Spline (spatial analyst) and Geostatistical Analyst?

texts:
Why the interpolation result of Spline with Tension in the Spatial Analyst is different from the that in the Geostatistical Analyst (both in ArcGIS 10), including their minimum values, maximum values, and isopach?

Just as the maps below, with the same original date and corresponding parameter setting, fig. 1 is done by Spline with Tension in Spatial Analyst, fig. 2 is done by Spline in Geostatistical Analyst, and they are totally different.

What does the prediction map in Geostatistical Analyst exactly mean?

As we know, Spline is an exact interpolation method, but in Geostatistical Analyst, it seems not. Just as the yellow dots  in the fig. 2, the original value of them are all 0, but why they appear in different color bands in the prediction map, which means that their predicted value are different with their original value.

Why? I'm really confused about that.

Some further explanations:

the yellow and red dots=0, and the values of the green dots are labeled in the figs, all of them are used to do the interpolation.
Fig. 1 is understandable to me, which conforms with definition of Spline, and all the dots appear in the correct color areas that include their original data.
I cann't understand fig. 2.

Why does the same interpolation method in different tools lead to very different results?

Forgive my poor English.

interpolation result of Spline with Tension in the Spatial Analyst

Prediction Map in the Geostatistical Analyst


texts: Using ArcPy to add polygon to FeatureClass?

texts:
I want to convert some poins to polygon and then  add the polygon to feature class.
but this does not work.where this code is wrong.

import arcpy
array = arcpy.Array([arcpy.Point(0, 0),
                  arcpy.Point(0, 1000),
                  arcpy.Point(1000, 1000),
                  arcpy.Point(1000, 0)
                  ])
polygon = arcpy.Polygon(array)
FC = arcpy.CreateFeatureclass_management("in_memory","Mylayer","POLYGON","","DISABLE   D","DISABLED")
with arcpy.da.InsertCursor(FC,("Shape@")) as iCur:
     iCur.insertRow(polygon)

Runtime error
Traceback (most recent call last):
File "<string>", line 2, in <module>
TypeError: argument must be sequence of values


texts:
The organisation I've just started at uses QGIS 2.14 Essen & I want to download it to practice at home on my Win 64 Home Premium machine. I've found an index of exe files (download.osgeo.org/qgis/win64/) but have no idea which one I should download.

e.g. what's the difference between this one QGIS-OSGeo4W-2.14.22-1-Setup-x86_64.exe & this one QGIS-OSGeo4W-2.14.22-1-Setup-x86_64.exe.md5sum? What's the md5sum part mean?


texts: How to measure the length/area of a line/polygon by just selecting it?

texts:
I can get what I want by using this code but have to do it everytime i select. The geometry of the features are small (city scale).

layer = qgis.utils.iface.activeLayer()
features = layer.selectedFeatures()
for f in features:
    geom = f.geometry()
    #print "length is ", geom.length()
    #print "area is ", geom.area()

Code taken from here -

Calculate line lengths with Python in QGIS


texts:
I'm tring to create a link to survey123 on browser.
I've been able to configure link for the app.
Link has become:
arcgis-survey123://?itemID=XXXXXXXX&field:ObjID=1&center=9.999,9.999

But web (Surveylink) version does not work. I can’t set center.
This does not work:
https://survey123.arcgis.com/share/XXXXXXX?field:ObjID=1&field:center=9.999,9.999

How can I set up my center for browser version?


texts: Where are the Spatialite settings stored?

texts:
after moving a spatialite-database to another directory, qspatialite says "Can't connect to Database" / Error unable to open database file. Even no project is open or I try to create a new project.
In Qgis browser it lists the spatialite-database but says "database does not exist".
I scanned all Qgis project files for the according database, but it does not show up.

The questions are:

where are these settings stored?
how to change these settings to the correct path?


texts: Using QGIS Rotate Point Symbols tool when multiple point objects have same coordinates

texts:
I'm using SVG symbols to display a layer which represents the location of traffic accidents.  The rotation of the symbol indicates the direction the vehicle was travelling.  The rotation angle is stored in a column of the table and the style uses this value to set the symbol rotation.  The layer also contains a column which indicates the general direction of travel, but it only contains one of four allowable values (North, South, East or West).  I have used this data to populate the symbol rotation column.
Now sometimes the angle needs some tweaking when the streets don't run exactly North-South.  I can edit the angle value directly, and sometimes I do that, but other times its convenient to use the Rotate Point Symbols tool:

My problem is when there are multiple accidents at the one location (typically an intersection) the  Rotate Point Symbols tool always rotates the top most symbol.  I can't find a way to drill down into the stack.
I have tried holding down various modifier keys.  I have tried selecting one of the lower objects via the attribute table, but the topmost still gets rotated.  I have style rule in the rendering order of the layer so selected objects are render on top - doesn't help.
The layer is from a third party and I would prefer not to modify the actual point objects so they don't overlap.
Looking for suggestions, thanks.


texts: How do I merge US counties consistently?

texts:
I am new to QGIS.

My project involves county-level data. The problem is that US counties keep changing! (See the US Historical Counties project: http://publications.newberry.org/ahcbp/)

Since my other data do not keep account of these changes, I plan to merge counties to create a single consistent geographic unit over time.

For example, if county A splits into county A and county B, then I simply treat county B as county A for the remainder of the time sample.

Are there any good ways to do this in QGIS, without arduously going through the attribute table?


texts: Reclassifying Rasters Error 001100 Failed because no statistics are available

texts:
I need to reclassify some vector data which was made from Euclidian distance from points and lines.

So after getting the Euclidian distances, I tried to run reclassify, but got:

  Error 001100 Failed because no statistics are available.

I ran BatchCalculate Statistics on all of the Rasters I want to reclassify, and that went through.

Tried Reclassify again, but failed, same error.

Went to try and use Get Raster Properties and got:

  Error 999999 Error Executing Function.

Here's my environments (I'm using ARCMAP 10.1)


texts: Creating a shapefile with SQL Query Using Geotools Hibernate and Java

texts:
Here is my context :

1] I have geometries in an postgresql database

2] My application uses JAVA/Hibernate

3] I Currently use Geotools to create a shapefile

How I do it ?

First i retrieve my results:

public List<Map<String, Object>> getObjectsToShape(String strQuery) throws BusinessException {

    Query query = entityManager.createNativeQuery(strQuery);

    org.hibernate.Query hbQuery = ((org.hibernate.jpa.HibernateQuery) query).getHibernateQuery();
    hbQuery.setResultTransformer(AliasToEntityMapResultTransformer.INSTANCE);
    List<Map<String, Object>> resultList = hbQuery.list();

    return resultList;
}

And example of a query :

    SELECT car.id AS ID,  \
       car.name as NAME, \
       car.creationdate as CR_DATE, \
       concat(car.id, car.name) as EXT_NAME, \
       car.notepercentage as NOTE_PER, \
       car.codepercentage as CODE_PER, \
       (select incidence from technicalmeasures where iddb=(select
               min(iddb) from technicalmeasures
               where technicalmeasures.car_iddb=car.iddb)) as INCID_ANG, \
       (select rotation from technicalmeasures where iddb=(select min(iddb) from technicalmeasures where technicalmeasures.car_iddb=car.iddb)) as ROTATION, \
       (select center from tiysd where iddb=(select min(iddb) from tiysd where tiysd.car_iddb=car.iddb))   as CENTER, \
       car.quality as QUALITY, \
       0.5 as MARK, \
       ST_AsText(car.shape) as shape  \
       from car

The query is more likely to be a view, many columns of other tables are selected, joinings, and static values..
The query itself is mandatory.
The map's key is the alias of the column name. and the value is the actual data.

the results contain different types :

The Geometry of course
Floats
Strings

Second,
I create my shapefile using "gt-shapefile" library of Geotools

// we only need one row to create fields of the shapefile
Map<String, Object> mapObject = objectList.get(0);

Set<String> keys = mapObject.keySet();

for (String columnName : columnNames) {
   try {
       wktReader.read(String.valueOf(mapObject.get(columnName)));
   } catch (ParseException exc) {
       builder.length(100).add(columnName.toUpperCase(),String.class); // first point i can not seem to retrieve the exact type because i don't know it, and can not find any solution without being too specific
   }
}

Then fill the shapefile :

      // Build feature
  SimpleFeatureType featureSchema = builder.buildFeatureType();
  // Feature Collection
  DefaultFeatureCollection featureCollection = new DefaultFeatureCollection();

  for (Map<String, Object> row : objectList) {

      SimpleFeatureBuilder simpleFeatureBuilder = new SimpleFeatureBuilder(featureSchema);

      for (String key : row.keySet()) {

          String value = String.valueOf(row.get(key)).trim();

          // Test if value starts with polygon on Multipolygon prefix
          if (value.startsWith("POLYGON") || value.startsWith("MULTIPOLYGON")) {
              wktReader.read(String.valueOf(value));
              simpleFeatureBuilder.set("the_geom", wktReader.read(String.valueOf(value)));
          } else {
              simpleFeatureBuilder.set(key.toUpperCase(), value);
          }

      }

      SimpleFeature simpleFeature = simpleFeatureBuilder.buildFeature(null);

      featureCollection.add(simpleFeature);

  }

My issues

This is not a satisfying method, I have searched for almost a year to reproduce the exact same work of "pgsql2shp". I looked deeper into Geotools libraries like  gt-jdbc.. but again could not find any working example. Is it possible to reproduce the exact work of pgsql2shp in Java ? with or Without Geotools ? if not, is it possible to not use hibernate and use only internal geotools librairies (a working example is very welcome)
In case of keeping hibernate. Do you know if it is possible to keep the column order from the sqlQuery ?


texts:
I need to filter out the low speed winds from the values inside a Grib file. This is my first time using this kind of data. I can read the file, and iterate inside all the GeoSpatialValues (getting lat, long and value from these).

To get the maximum resolution, should I use any specific list of the so called "component of wind" ?

Is there an easy way to calculate the speed from these values or should I iterate in the U values and then V values of wind to finally combine them into a vector?


texts: How can I display an image on QGIS from a DJI drone using python?

texts:
I have a few images from a DJI drone and I can get all of the information from the image exifread.   I am trying to get an image to be displayed in the correct location on a QGIS map.  Is there an example of how to use this data to get the image to show correctly on a QGIS map?

I have tried the GDAL and nothing happens yet.

https://www.gdal.org/gdal_tutorial.html

Elements
drone-dji:AbsoluteAltitude= 258.91
drone-dji:RelativeAltitude= 60.9
drone-dji:GimbalRollDegree= 0.0
drone-dji:GimbalYawDegree= -148.3
drone-dji:GimbalPitchDegree= -90.0
drone-dji:FlightRollDegree= -2.8
drone-dji:FlightYawDegree= 179.4
drone-dji:FlightPitchDegree= -7.3

Exif Data
Key: Image ImageDescription  value  DCIM\100MEDIA\DJI_0002.JPG
Key: Image Make  value  DJI
Key: Image Model  value  FC550
Key: Image Orientation  value  Horizontal (normal)
Key: Image XResolution  value  72
Key: Image YResolution  value  72
Key: Image ResolutionUnit  value  Pixels/Inch
Key: Image Software  value  v01.29.5379
Key: Image DateTime  value  2019:04:25 10:50:33
Key: Image YCbCrPositioning  value  Centered
Key: Image ExifOffset  value  182
Key: GPS GPSVersionID  value  [0, 0, 0, 0]
Key: GPS GPSLatitudeRef  value  N
Key: GPS GPSLatitude  value  [42, 31, 524857/10000]
Key: GPS GPSLongitudeRef  value  W
Key: GPS GPSLongitude  value  [83, 6, 552361/10000]
Key: GPS GPSAltitudeRef  value  0
Key: GPS GPSAltitude  value  258907/1000
Key: Image GPSInfo  value  710
Key: Image XPComment  value  [48, 0, 46, 0, 57, 0, 46, 0, 49, 0, 57, 0, 57, 0, 0, 0, 0, 0, 0, 0, ... ]
Key: Image XPKeywords  value  [78, 0, 45, 0, 65, 0, 83, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Key: Thumbnail Compression  value  JPEG (old-style)
Key: Thumbnail XResolution  value  72
Key: Thumbnail YResolution  value  72
Key: Thumbnail ResolutionUnit  value  Pixels/Inch
Key: Thumbnail JPEGInterchangeFormat  value  41972
Key: Thumbnail JPEGInterchangeFormatLength  value  8148
Key: EXIF ExposureTime  value  521/500000
Key: EXIF FNumber  value  11/5
Key: EXIF ExposureProgram  value  Program Normal
Key: EXIF ISOSpeedRatings  value  110
Key: EXIF ExifVersion  value  0230
Key: EXIF DateTimeOriginal  value  2019:04:25 10:50:33
Key: EXIF DateTimeDigitized  value  2019:04:25 10:50:33
Key: EXIF ComponentsConfiguration  value  CrCbY
Key: EXIF CompressedBitsPerPixel  value  6343987/1990656
Key: EXIF ShutterSpeedValue  value  -9906/-1000
Key: EXIF ApertureValue  value  227/100
Key: EXIF ExposureBiasValue  value  0
Key: EXIF MaxApertureValue  value  17/10
Key: EXIF SubjectDistance  value  0
Key: EXIF MeteringMode  value  Spot
Key: EXIF LightSource  value  Unknown
Key: EXIF Flash  value  No flash function
Key: EXIF FocalLength  value  15
Key: EXIF FlashPixVersion  value  0010
Key: EXIF ColorSpace  value  sRGB
Key: EXIF ExifImageWidth  value  4608
Key: EXIF ExifImageLength  value  3456
Key: Interoperability InteroperabilityIndex  value  R98
Key: Interoperability InteroperabilityVersion  value  [48, 49, 48, 48]
Key: EXIF InteroperabilityOffset  value  680
Key: EXIF ExposureIndex  value  0/0
Key: EXIF FileSource  value  Digital Camera
Key: EXIF SceneType  value  0
Key: EXIF CustomRendered  value  Normal
Key: EXIF ExposureMode  value  Auto Exposure
Key: EXIF WhiteBalance  value  Auto
Key: EXIF DigitalZoomRatio  value  0/0
Key: EXIF FocalLengthIn35mmFilm  value  30
Key: EXIF SceneCaptureType  value  Standard
Key: EXIF GainControl  value  None
Key: EXIF Contrast  value  Normal
Key: EXIF Saturation  value  Normal
Key: EXIF Sharpness  value  Normal
Key: EXIF DeviceSettingDescription  value  [0, 0, 0, 0]
Key: EXIF SubjectDistanceRange  value  0
Key: EXIF LensMake  value  UnKnown
Key: EXIF LensModel  value  DJI MFT 15mm F1.7 ASPH
PS C:\Users\Cary\Documents\Programming\Python\RailProfile>


texts: Check for mxd write access with python addin

texts:
I'm looking for a way to check for write access to an mxd upon the initiation of a python addin tool. My tool should only be enabled if the user has write access to the mxd. I thought there was a way to check for read/write access using arcpy (or maybe os?) but it's lost on me now. I tried the following as a work-around (saving upon initiation):

class ToolClass2(object):
    def __init__(self):
        mxd = arcpy.mapping.MapDocument ("CURRENT")
        try: mxd.save ()
        except IOError:
            pythonaddins.MessageBox ("The map document is locked.",
                                     "Map file locked", 0)
            self.enabled = False

However, this doesn't produce an IOError. Instead I get a warning box (below) from ArcMap and the try is successful and the except is ignored.

I've also tried making use of os.access, but without results (os.access (mapPath, os.W_OK) returns True):

mapPath = mxd.filePath
if not os.access (mapPath, os.W_OK):
    pythonaddins.MessageBox ("The map document is locked.",
                             "Map file locked", 0)
    self.enabled = False

I've now realized that I was thinking of arcpy.TestSchemaLock, but according to the help this is only for feature classes, tables, and feature datasets.

Developing with version 10.4 but the end users will be 10.1.


texts: Georeferencing Moon photos and LRO DEM but how to model the curvature of the Moon?

texts:
I tried to search for solution to this, but I could not find it. However, if this matter has been dealt somewhere here, please forward me to a correct thread.
To summarize:
I am using QGIS and Qgis2threejs plugin to create 3D models from mine and my astroimager colleagues' photos of the Moon utilizing LRO DEM. However, I lack the curvature of the Moon in my models and I would like to add it if possible.
Software:
QGIS 3.10.13 (Coruna, LTR), Qgis2threejs 2.5
Details:
I utilize LRO mission DEM (https://astrogeology.usgs.gov/search/details/Moon/LRO/LOLA/Lunar_LRO_LOLA_Global_LDEM_118m_Mar2014/cub) and Lunar photos taken e.g by me. This way I can easily do georeferencing between my images and LRO DEM and create a 3D models running in web browser showing Lunar surface features and elevation. However, with large Moon mosaics it would be nice to get the curvature of the Moon taken into account so that the horizon would curve instead of the horizon being completely flat. I do not know how / if Moon curvature can be parameterized into the model.
Examples:
(1) This shows a model in which almost a complete strip of the Moon along the terminator from south to north is modeled in 3D using QGIS and Qgis2threejs. The model shows the surface features (elevation) but lacks the curvature of the Moon which would be clearly seen in this kind of a large area model. Please note that the z-axis (elevation) has been scaled with a factor of 1.5 to show the elevation in a bit more visual way. The model may take some time to load (typically done within 1 minute or less) but should be working with many browsers. Model can be opened here.
(2) This is provided as a reference only: A smaller area model does not require modeling the curvature of the Moon as these kind of models are OK as they are. However, if possible, also this model could possibly benefit if the curvature could be modeled. This model can be opened here.
What have been tried:
I have no idea how to implement the Moons's curvature into my models. I have searched a lot but I haven't found a solution for the issue. Thus, I have to ask from those who are more capable users than me. Please note that I only know the basics of QGIS and am not in any way an advanced user. I read e.g. this, but I am not sure how to e.g. modify the DEM (LRO DEM geotiff) to take into account the Moon's curvature.
What am I looking:
A solution, hopefully within QGIS or Qgis2threejs, to take into account the curvature of the moon into my 3D models generated via Qgis2threejs and exported to be run via web browser. I am not looking a solution that requires e.g. Blender or similar to create the curvature of the moon unless I still can utilize the web browser based representation which Qgis2threejs exports (which are pretty awesome). Hopefully such a process of taking the curvature of the Moon into account could be accomplished within QGIS.


texts:
I'm trying to merge 12 .tif raster layers (raster -> misc -> build virtual raster), but keep getting this error:

  Loading resulting layers The following layers were not correctly
  generated./var/folders/4f/3vtn5c5s7zd3qnrn03_5bc1r0000gn/T/processing_201ee1195d8f460f8c74a248f1bab91c/878b8a3be722499da52c93b6f72c083c/OUTPUT.vrtYou
  can check the 'Log Messages Panel' in QGIS main window to find more
  information about the execution of the algorithm."

  "2018-10-19T13:21:01     INFO    gdalbuildvrt -resolution average
  -separate -r nearest -input_file_list /var/folders/4f/3vtn5c5s7zd3qnrn03_5bc1r0000gn/T/processing_201ee1195d8f460f8c74a248f1bab91c/buildvrtInputFiles.txt
  /var/folders/4f/3vtn5c5s7zd3qnrn03_5bc1r0000gn/T/processing_201ee1195d8f460f8c74a248f1bab91c/523821c6a9a64cd2a7789f266b312785/OUTPUT.vrt
  2018-10-19T13:21:01     INFO    GDAL execution console output
               /bin/sh: gdalbuildvrt: command not found


texts: How to Select By Location with matching attribute?

texts:
I have a whole bunch of overlapping polygons (buffers), I want to select with points only those poly's which:

the point falls within; and
the poly attribute matches the point attribute

What would be the best way to do this selection, perhaps with ModelBuilder? (using ArcGIS Desktop 9.3.1)


texts: Drawing Manager circle moves when zoom out

texts:
I am using DrawingManager to draw a circle at zoom level 19.

this.drawingManager.setDrawingMode(google.maps.drawing.OverlayType.CIRCLE);

when zoomed out the circle is moving. Now the circle is in the intersection.

Is there any way to prevent this?


texts: I would like to know if there are any open source flex viewers some thing on the same line as ArcGIS Flex Viewer  ??


texts: Seeking marketplace intended to commercialize QGIS applications?

texts:
If we want to sell general application such as website theme, plugins, mobile app, game, etc, we can use envato marketplace, google play store, or other marketplaces for selling and buying application.

However, i need to knows if there is a marketplace that functioning to sell and buy application in the realm of GIS such as Webgis template, QGIS plugin, or standalone application used for geospatial needs, etc?


texts: I use OpenLayers and GeoServer. I can read kml file format from OpenLayers. And I want to edit this polygon in kml. So I want polygon data(polygon circle and center coordinate etc..) to GeoServer.


texts:
I'm having some trouble figuring out the best way handle some GIS data. I want to analyze pollution from specific sites at specific points. I have a list of sites I'm interested in all across the continental U.S. And I have a method for getting pollution data for those sites. The pollution data is defined by the distance from the site, so I've decided that I need to store that data in a polar grid. I thought it would be interesting to figure out how to find the average amount of pollution for a given point. I need to find the average amounts of pollution for points that that are in multiple intersecting polar grids.

I have:

a few sites across the continental U.S. in the same projection.
360x360 polar grids with a radius of 200m. (I'm not attached to those numbers, but the grid size will be constant across all grids and the radius needs to be exact and somewhere around 200m)
polar grids with the average size of about a square mile.

And I want to analyze pollution at points in intersecting polar grids.

My last attempt was to use an Albers Equal Area Conic projection and then create polar grids from each point and find the resulting geometric shapes from the intersecting polar grids. Then for any given point in a stitched together area I could find out what pollution sources are within 200m.

But I'm trying to get a handle on how accurate this will be. First my test polar grids are not exactly circular with a radius of 200m. Distances from center to edge are usually within 10%, but they aren't exact. This is a bigger deal if I'm wanting to deal with intersecting polar grids of size 360x360. This means that if distance is off drastically I will get errors in my intersecting polygons. But I also want to map specific points on the ground to intersecting polygons so I want to have an accurate area. The other problem is the size, I briefly looked at individual state plane zones, but I would be using multiple zones if I was stitching together multiple pollution sources.

Any suggestions on how to handle this? It's looking like a really messy problem. My main goal is to take be able to take any point in the coverage area and get all polar grid cells that contain that point.


texts: Sorting through a column in Postgres which is not numeric

texts:
I am using postgres and i want to sort the data using order by but there is an issue with the sorting which i am unable to figure out.

I have a column which is character varying and the values in that column are.

_________
|ID  |Name|
|1   | ABC|
|5   | ABC|
|2   | ABC|
|7-A | ABC|
|2-A | ABC|
|3   | ABC|
-----------

I want the results as following

_________
|ID  |Name|
|1   | ABC|
|2   | ABC|
|2-A | ABC|
|3   | ABC|
|5   | ABC|
|7   | ABC|
-----------


texts: Export vector layer to sqlite database using PyQGIS

texts:
I have successfully created a (virtual) polygon layer named vlyr. I also have a QgsFeature object called poly. I also added two fields to my layer called "name" and "size".

Now I want to export my polygons to a sqlite database which I manually created for testing purposes. It is called test.sqlite.

I am able to establish a db connection and create a new sqlite table using the following code

import sqlite3 as lite
import sys

con = lite.connect('test.sqlite')
for l in vlyr.getFeatures():

    with con:
        cur = con.cursor()
        abc = (
              (1, 'fstName', 6543.12),
              (2, 'sndName', 12635.11)
            )
        cur.execute("DROP TABLE IF EXISTS Polygonx")
        cur.execute("CREATE TABLE Polygonx (Id PRIMARY KEY, Name TEXT, Size REAL)")
        cur.executemany("INSERT INTO Polygonx VALUES (?,?,?)", abc)
        cur.execute("SELECT * FROM Polygonx")
        rows = cur.fetchall()
        for row in rows:
            print row

I also stored the geometry in a variable called a, but am currently not doing much with it:

a = (l.geometry().exportToWkt())

My question is: how can I create abc dynamically, that is, from my existing polygon layer and add the geometry?


texts:
I'm having a problem when setting the IactivewViews extent. I have done it the same way on other places and it does work. In here, it just throws me a LoaderLockException.

        // OLD CODE
        private void airspaceElementsTreeView_NodeMouseDoubleClick(object sender, TreeNodeMouseClickEventArgs e)
        {
            if (e.Node.Name.StartsWith("feature_"))
            {
                IFeatureWorkspace workspace = ((IDataset)_LoadedAirspace.UnderlyingObject.Table).Workspace as IFeatureWorkspace;
                ITable relatedTable = workspace.OpenTable(e.Node.Parent.Name);

                int oid = GetNodeObjectId(e.Node);
                IRow tempRow = relatedTable.GetRow(oid);

                IGeometry geom = null;

                if (tempRow is IFeature)
                {
                    IFeature feat = tempRow as IFeature;

                    geom = feat.ShapeCopy;
                    ISpatialReference srCopy = ((IClone)geom.SpatialReference).Clone() as ISpatialReference;

                    IEnvelope envelope = new EnvelopeClass();
                    envelope.SpatialReference = srCopy;
                    envelope.Union(geom.Envelope);
                    envelope.Expand(.25, .25, false);

                    // error here!
                    IActiveView view = ((IMxDocument)m_application.Document).ActiveView;
                    view.Extent = envelope;
                }
                else
                    return;

                Marshal.ReleaseComObject(relatedTable);
            }
            else
                return;
        }

By the way, this is inside a IDockableWindowDef. Any chance that THAT is ruining things for me?

edit following the suggestions of my fellow friends, I did on two separate statements. First I'm trying to get IMxDocument from m_application and then casting it to IActiveView.

It crashes when I'm trying to cast m_application.Document to IMxDocument (check new code below)

    private void airspaceElementsTreeView_NodeMouseDoubleClick(object sender, TreeNodeMouseClickEventArgs e)
    {
        if (e.Node.Name.StartsWith("feature_"))
        {
            IFeatureWorkspace workspace = ((IDataset)_LoadedAirspace.UnderlyingObject.Table).Workspace as IFeatureWorkspace;
            ITable relatedTable = workspace.OpenTable(e.Node.Parent.Name);

            int oid = GetNodeObjectId(e.Node);
            IRow tempRow = relatedTable.GetRow(oid);

            IGeometry geom = null;

            if (tempRow is IFeature)
            {
                IFeature feat = tempRow as IFeature;

                geom = feat.ShapeCopy;
                ISpatialReference srCopy = ((IClone)geom.SpatialReference).Clone() as ISpatialReference;

                IEnvelope envelope = new EnvelopeClass();
                envelope.SpatialReference = srCopy;
                envelope.Union(geom.Envelope);
                envelope.Expand(.25, .25, false);

                // crash here
                IMxDocument document = m_application.Document as IMxDocument;
                IActiveView view = document.ActiveView;
                view.Extent = envelope;
            }
            else
                return;

            Marshal.ReleaseComObject(relatedTable);
        }
        else
            return;
    }


texts:
I've just updated to QGIS Brighton on Ubuntu 14.04, and I've noticed that when I open Openlayers plugin, in order to load google layers, the whole program suddenly crashes.

The only log I see, when open with terminal is this:

setLayerType: Google Hybrid
 setLayerType: Google Hybrid
OpenlayersController.__init__
[WORKER THREAD] Calling request() asynchronously
[WORKER THREAD] Waiting for the async request to complete
[GUI THREAD] Processing request
page file: file:////home/pablo/.qgis2/python/plugins/openlayers_plugin/weblayers/html/google_hybrid.html
OpenlayersWorker waiting for page to load
Violación de segmento (`core' generado)

Does anybody know something about?


texts:
at the moment I am trying to create an indexed vector grid via "vector -> research tools -> vector grid...", in the following format:

A1 A2 A3 ...
B1 B2 B3 ...
C1 C2 C3 ...
.. .. ..

The creation of the grid is no Problem, but I am stuck with labeling the individual polygons.

I already wrote an Excel Visual Basic Project, which creates a script for the field calculator and helps me to label each row from A to X. Now I want to label the columns with numbers from the first one to the last. Do you have any idea how I could use the @row_number to specify the column?

Here is the output of the script for the field calculator:

CASE
WHEN ID>-1 AND ID<5 THEN concat('A',@row_number)
ELSE
CASE
WHEN ID>4 AND ID<10 THEN concat('B',@row_number)
ELSE
CASE
WHEN ID>9 AND ID<15 THEN concat('C',@row_number)
ELSE
CASE
WHEN ID>14 AND ID<20 THEN concat('D',@row_number)
ELSE
CASE
WHEN ID>19 AND ID<25 THEN concat('E',@row_number)
END
END
END
END
END


texts: Can you search WMS services for a particular projection?

texts: Is there an interface or website to search public WMS services for location/projection?


texts: Finding Band Count in Rasters using ArcPy?

texts:
I'm having troubles finding the band count for my raster files using ArcPy.

I need to first list all the raster files in a directory, and then display the band count for each file. I am able to list the files, but when it gets to the band count, it gives me an error stating that Method Band Count does not exist.

The following is what I have for my code so far:

import os
import arcpy

arcpy.env.workspace = "C:\\Rasters" # Contains TIFF, IMG, GRID formats
filePath = arcpy.env.workspace

rasterList = arcpy.ListRasters("*", "ALL")
desc = arcpy.Describe(filePath)

for names in rasterList:
    print names
    print desc.bandCount

Any ideas where I'm going wrong? I was using the Geoprocessor Programming Model as a reference.


texts: Create a path in mapping software and export as vector?

texts:
Is it possible to map a route in any mapping software (OpenStreetMaps, Google Maps/Earth, QGIS, etc. preferably Open Source or free) and export this route as a vector (PDF, AI, SVG etc.)?

I'm only interested in the route taken (as a line), without any other information.


texts: Merging multiple feature classes into one feature class

texts:
I have several feature classes in different geodatabases and I am trying to merge into one feature class, same geometry one attribute table in one geodatabase (attribute fields are the same).
It keeps coming out as a shapefile which I cannot have as it gets rid of my domains. How do I merge and have a resulting single feature class and/or is there something else I should use other than the merge tool?
Background: Several feature classes originally in their own geodatabase. I copied them to a single geodatabase. The feature classes are polyline culvert data originally collected with a Trimble GPS Unit.
I do not know how to script or use Python yet so I am trying to do this using the GUI.


texts:
How to set Layer's options' title? I want to identify each layer

 var overLayers = [
 {
   name: "Bar",
   icon: iconByName('bar'),
   layer: L.geoJson(Bar,{onEachFeature: onEachFeature}),
   group: "Bar",
   marekerType:"accordion",
   active:true,
 }
];


texts: Getting error after trying to harvest from geoserver into pycsw

texts:
I'm using osgeo live to test the pycsw functionality. I need to harvest metadata from an existing Geoserver service (previously configure for this task).

I'm quite new with CSW and I'm not sure if my procedure is right. I've trie to follow the instrutions under:
http://geopython.github.io/pycsw-workshop/docs/intro/intro-exercises.html#metadata-harvesting

My xml config file (called 'cat.xml') for the harvest is:

  <?xml version="1.0" encoding="UTF-8"?>
  <Harvest xmlns="http://www.opengis.net/cat/csw/2.0.2" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/cat/csw/2.0.2 http://schemas.opengis.net/csw/2.0.2/CSW-publication.xsd" service="CSW" version="2.0.2">
    <Source>http://localhost:8082/geoserver/wms</Source>
    <ResourceType>http://www.opengis.net/wms</ResourceType>
    <ResourceFormat>application/xml</ResourceFormat>
  </Harvest>

for the source element I've used this one too (without exit):

<Source>http://localhost:8082/geoserver/cat/wms</Source>

I use the command:

$ cd /var/www/pycsw
$ python ./sbin/pycsw-admin.py -c post_xml -u http://localhost/pycsw/csw.py -x ~/cat.xml

Then I get the following error message:

Intializing static context
Executing HTTP POST request /home/user/cat.xml on server http://localhost/pycsw/csw.py
Traceback (most recent call last):
  File "./sbin/pycsw-admin.py", line 245, in <module>
    print admin.post_xml(CSW_URL, XML, TIMEOUT)
  File "/usr/lib/python2.7/dist-packages/pycsw/admin.py", line 506, in post_xml
    raise RuntimeError(err)
RuntimeError: HTTP Error 401: Unauthorized

The request GetCapabilities works:

localhost/pycsw/csw.py?config=tests/suites/default/default.cfg&service=CSW&version=2.0.2&request=GetCapabilites

Can you help me?


texts: How to store and serve the time a GeoTIFF refers to

texts:
I'm writing software that creates temperature maps as GeoTIFF files. I need to store the time to which the GeoTIFF file refers; not the raster creation time, but the time during which the temperature was as the raster data says.

What I did was store it in the GeoTIFF file as a GDAL metadata item; i.e. in a GDALMetadata private tag. The metadata item is called "TIMESTAMP".

I now need to serve the GeoTIFF file as a WCS with MapServer.

Is storing the time as a GDAL metadata item a good idea? Is it standard practice?
How is this going to be accessed with WCS? Through DescribeCoverage?
How can MapServer be told to get the metadata item from the GeoTIFF and somehow serve it to the client?


texts: Debugging HTML/JavaScript code for OpenLayers?

texts:
I'm really stumped on why my HTML/Javascript code isn't working. My webpage currently only shows my map div, but no actual map content.

My code so far:

HTML:

Javascript:
var map;

function initMap() {
 //initialise variables
 var osmlayer;
 var google_terrain;
 var point;

 //create new map object
 map = new OpenLayers.Map("mymap");

 //create new API layers
 osmlayer = new OpenLayers.Layer.OSM("Open Street Map");
 google_terrain = new OpenLayers.Layer.Google("Google Terrain Map", {
     type: google.maps.MapTypeId.TERRAIN,
     numZoomLevels; 22
 });

//add layers to the map
 map.addLayer([osmlayer, google_terrain]);

 map.addControl(new OpenLayers.Control.LayerSwitcher());

 // define point as a new LonLat object and transfom
 point = new OpenLayers.LonLat(-6.67, 55.15);
 point.transform("EPSG:4326", "EPSG:900913");

 //center the map and set zoom level
 map.setCenter(point,10);

 }


texts: GRASS scripting for producing multi-band raster?

texts:
I am trying to loop the r.sun function within a Python script.

Anyway, I seem not to understand how I can stack the raster maps each loop would produce into a single multi-band raster (preferably in geotiff format).

I think one way would be using the list produced with g.list or (even better) g.mlist as input parameter for the i.group function and then export the created group as multi-band file (as shown in the example here).
Here is my simple code, modified from this page.

# import the needed modules
from osgeo import gdal
import numpy as np
import grass.script as grass

# prepare the needed input
dem = "dem4comparison@globrad_test"
doy = np.arange(1,366) # day of the year

# calculate slope and aspect maps
grass.run_command("r.slope.aspect",
              elevation = dem,
              slope = "dem.slope@globrad_test",
              aspect = "dem.aspect@globrad_test")

# i is the day of year (DOY) to use for running r.sun
for d in doy:
    # want to get global radiation
    radname = "globrad" + str(d)
    # run the solar model
    grass.run_command("r.sun", flags = "s",
              elevin = dem,
              aspin = "dem.aspect@globrad_test",
              slopein = "dem.slope@globrad_test",
              insol_time = insolname, glob_rad = radname, day = d, step = 1)
    # export the outputs to GTiff files
    grass.run_command("r.out.gdal", input = radname + "@globrad_test", output = radname)

The last row would export each map to a geotiff file after each loop. Rather, I would stack them all into the same multi-band file with n*m*p dimensions, where m and n are rows and columns of my images (each of the same size of course) and p the number of output maps after the loops (that is, 365, the number of the days in a year).


texts: Is there a way to prevent labels from 2 different layers (province and city names)? This map of Argentina looks a bit messy as some labels are overlapping and become unreadable.


texts:
How do you convert a JSON file to a GeoJSON file or csv. I have a JSON file with lat lon information, but how can I can convert this file to GeoJSON or csv. I have found some online tools to convert it to csv, but I want to use it in a batchfile, so I need a tool.

The JSON file looks like this:

[
{"date":"2014-09-25","time":"20:49:09","lat":"53.269","lon":"6.935","depth":"3.0","place":"Meedhuizen","mag":"1.5","type":"GH","last":"True"},
{"date":"2014-09-24","time":"23:49:36","lat":"53.351","lon":"6.693","depth":"3.0","place":"Huizinge","mag":"1.0","type":"GH","last":"False"},
{"date":"2014-09-23","time":"17:34:48","lat":"53.242","lon":"6.728","depth":"3.0","place":"Lageland","mag":"1.0","type":"GH","last":"False"}
]

Hope someone can help


texts:
[QGIS V.3.10.3 A Coruña]

For streets in my fantasy world map I want to use a line ESRI shapefile.

My custom CRS is a lonlat one and looks like the following:

GEOGCRS["unknown",
  DATUM["unknown",
    ELLIPSOID["Yartol",6370997,326.249334289226,
      LENGTHUNIT["metre",1,ID["EPSG",9001]]]],
PRIMEM["Kangorr",0,
  ANGLEUNIT["degree",0.0174532925199433,
    ID["EPSG",9122]]],
CS[ellipsoidal,2],
AXIS["longitude",east,
  ORDER[1],
  ANGLEUNIT["degree",0.0174532925199433,
    ID["EPSG",9122]]],
AXIS["latitude",north,
  ORDER[2],
  ANGLEUNIT["degree",0.0174532925199433,
    ID["EPSG",9122]]]]

For my understanding LENGTHUNIT["metre"... should do the stuff. The measuring bar and tools show good results. In my project properties everything is set to give me correct measurements when scrolling around and measuring stuff, but the lines of my shapefile fill up the whole screen when I set the line thickness to 1 meter at scale, even when I should see my whole world map! The CRS is set both for the project and the shapefile layer. The (apparently unchangeable) properties of the layer indicate, that the unis is degrees though.

What am I missing here?


texts: How to load tiles of a specific area in OpenLayers 3 using a bounding box

texts:
I have built a simple OpenLayers 3 map, see below for code:

HTML code:

CSS:

#map{
     width:100%;
     height:800px;
}

OpenLayers map initialization code:

var map;
function init(){
    map = new ol.Map({
        target:'map',
        renderer:'canvas',
        view: new ol.View({
            projection: 'EPSG:900913',
            center:[-8015003.33712,4160979.44405],
            zoom:5
        })
    });
    var newLayer = new ol.layer.Tile({
    source: new ol.source.OSM()
    });
    map.addLayer(newLayer);
}

Now my map loads just fine, what I would really like to do is, I would like to have only a specific area loaded, for which I need a bounding box, I would like to load the the bounding box for THIS area. But I am new to OpenLayers 3, so I tried and I just couldn't make the bounding box.

Can somebody help me make the bounding box please?


texts: Saving rasters in a Python for loop fails only on last iteration

texts:
I'm a newcomer to Python scripting and I'm having what's seemingly a very simple problem but one that's stumped me enough to send me posting to these forums for the first time.

I'm using a for loop to save rasters to a specified folder. The loop executes perfectly, but the last iteration of the loop invariably produces an empty raster with no spatial reference. So, iteration 1 will be an empty raster until I run iteration 2, at which point the product of iteration 1 mysteriously transforms into what I want, & so on w/ iterations 2 & 3, etc. In keeping w/ this pattern, the last raster created by the loop remains empty and unprojected.

Code is as follows ("costsurf" & "distDirString" are defined in an earlier part of the script):

for Num in Nums:
 focal = arcpy.sa.SetNull("Corr","Corr","VALUE <> " + str(Num))
 costDist = arcpy.sa.CostDistance(focal,costsurf)
 withhole = arcpy.sa.SetNull(~ arcpy.sa.IsNull(focal),costDist)
 withhole.save(os.path.join(distDirString,"dist_%s.img" %Num))

It's supposed to pull out a chunk of the raster "Corr" with cell values Num by setting all else to NoData, run the Cost Distance tool on it, then poke a hole in the result where the source cells are. I have no idea how to think about what's failing here, since the loop does exactly what I want it to, but only if it's not the last iteration to be run.

It seems to behave itself if I get rid of the ".img" and save them as GRIDs -- but only sometimes?? And I'd really rather stick to .img if at all possible, which it surely must be. I'm using ArcGIS 10.0 btw. Thanks for your help!


texts: Buffer around shoreline hitting user memory limits

texts:
I am having issues hitting user memory limits creating a buffer around a shoreline. I think its because there are too many elements in the feature collection. This might be a tiling issue? How do I create a polygon or set of polygons that represent this buffer around this shoreline feature? Here is the code and asset shared:
var shoreline = ee.FeatureCollection("users/nateme16/eeshore");
Map.addLayer(shoreline)

var shoreline2=shoreline.geometry()
var buffer= shoreline2.buffer(200)

Map.addLayer(buffer)

link code


texts: Problems with 'layer exists ' while running ogr2ogr

texts:
I run command:

ogr2ogr -f PostgreSQL PG:"dbname=regional host=localhost user=*** password=***" Hsp_09_13.MIF -t_srs EPSG:3857 -nlt MULTIPOLYGON -overwrite -nln postcode_sector_boundaries

And I have error:

  ERROR 1: Layer postcode_sector_boundaries already exists, CreateLayer
  failed. Use the layer creation option OVERWRITE=YES to replace it.

So I run this command then:

ogr2ogr -f PostgreSQL PG:"dbname=regional host=localhost user=*** password=***" Hsp_09_13.MIF -t_srs EPSG:3857 -nlt MULTIPOLYGON -lco OVERWRITE=YES -overwrite -nln postcode_sector_boundaries

and got error:

  ERROR 1: CREATE TABLE "public"."postcode_sector_boundaries" ( OGC_FID
  SERIAL, PRIMARY KEY (OGC_FID) ) ERROR:  relation
  "postcode_sector_boundaries" already exists

So I remove table postcode_sector_boundaries and run once more last command and got:

  ERROR 1: 'CREATE INDEX "postcode_sector_boundaries_geom_idx" ON
  "public"."postcode_sector_boundaries" USING GIST ("wkb_geometry")'
  failed for layer postcode_sector_boundaries, index creation has
  failed.

postcode_sector_boundaries_geom_idx is unique key so I think that maybe my data are broken. I would like that ogr2ogr create table without that index but I can not find where definition for this create query can be. Any one knows? Or maybe there is another solution for my problem?

edit:

Ok, I know what was the problem - it was because there were some other tables using values from one I wanted recreate.


texts:
I have this raster wich represent some areas of observation from 3 given points following a series of rules considering elevation and direction of the slope.
Basically, the white-ish color represent no point can see this area (0); the orange area is seen from 1 point (1); and the brown area is seen by 2 points (2).
Now, I have to calculate the area of each color, specialy the orange and brown ones, but I don't have a clue on how to do it.
I tried to convert this raster to vector to calculate te area of the polygon; but when I do it there're lots of polygons and none of them seem to conservate the value of the pixel.

The raster I'm interested in was obtained by the cumulative viewshed tool.

Any idea on how to do it?


texts: Using Workflow Manager to log task times

texts:
I'm currently updating our Workflow Manager workflow used in production.  I'd like to add a way for supervisors to track hours used on different tasks.  For instance, after an edit task is completed, the users will be prompted to enter time spent on editing.  Similarly, after a QC task is complete, the user is prompted to enter time spent on QC.  Then at the end of the day, supervisors can open a query that summarizes hours spent on each task for each job/job group/etc.

I've added a few questionnaire steps to get users to manually input their time for specific tasks.  I am aware that the results of these are stored in the JTX_HISTORY table, however, the results are stored in a manner which is (to my eye) not easily queryable for a report in WM Administrator.  Is there a method to achieve what I am trying to accomplish?


texts: Best way to visualize the PostGIS Raster in OpenLayers

texts:
I have been stuck by this problem for a long time. Several ways have already been tried, but none of them works as I was expecting. Here are the description of each of these method as well as the problem I am facing.

I want to know which is the best option to follow on because I am getting confused now.

Thanks.

Method 1:  (GDAL + OpenLayers.Layer.Image);

Convert the PostGIS raster to a tif file first using gdal_translate, and then translate this tif file to image (PNG) file. (Because it is faster than convert it directly from PostGIS to PNG).

Then, load this image to Openlayers.
Problems: Don't know how to make it colorful. Because the default iamge is in gray color. Don't know how to add a color table to this.
Method 2:  (MapServer + OpenLayers.Layer.WMS):
Already get it work in a visual server with 2GB ram. But the problem is the speed. It takes long time to load the data on the map compared with load it as a image file.

Is that the problem with our virtual server or how the way we set up the MapServer?
Method 3: GeoServer + Openlayers.Layer.WMS:

Cannot get it work yet. Because I don't know how to connect it with PostGIS raster. On the previous Post, it seems we have to connect it to JDBC. But I didn't know someone has ever made it work or not.

Waiting for the help. Thanks.


texts: Starting GeoServer on OSGeoLive automatically

texts:
I want to use OSGeo Live VM 8.0 as a test GeoServer. Following its website, I started the server by

  Select the Geospatial ‣ Web Services ‣ GeoServer ‣ Start GeoServer in
  the menu

But after rebooting the server, geoserver is down again.

How do I start geoserver automatically as if I installed geoserver myself on a fresh Ubuntu?


texts: Buffering a line with varying distances associated with directional quadrant by QGIS?

texts:
I would like to buffer a set of lines with different distances to the NW, SW, SE, and NE. I have the distances stored in four attribute fields that correspond to the direction the buffer should be applied in. I have reviewed the two posts How to buffer a line with different distances for left/right sides and QGIS buffer polygon with different widths on each side but they do not fully address my issue. My lines do not flow East to West or North to South so the right/left side buffer option does not offer a solution.
Edit:I have attached a rough rendering of what I would like the output to be. The edges should be smooth.


texts:
I'm trying to find the closest point on a line to another point as in this post. However when I try to write the result as a shapefile I get an error:

  TypeError: QgsVectorFileWriter.addFeature(QgsFeature,
  QgsFeatureRendererV2 renderer=None, QGis.UnitType
  outputUnit=QGis.Meters): argument 1 has unexpected type 'Point'

import shapely
from shapely.wkt import *
from shapely.geometry import *
from qgis.gui import *
from qgis.core import *
from PyQt4.QtCore import Qt

lineLayer = iface.mapCanvas().layer(0)
pointLayer =  iface.mapCanvas().layer(1)
canvas =  iface.mapCanvas()
spIndex = QgsSpatialIndex()
lineIter =  lineLayer.getFeatures()

provider = pointLayer.dataProvider()
fields = provider.fields()
writer = QgsVectorFileWriter('Test.shp', 'CP1250', fields, provider.geometryType(), provider.crs(), "ESRI Shapefile")
for lineFeature in lineIter:
    spIndex.insertFeature(lineFeature)
pointIter =  pointLayer.getFeatures()
for feature in pointIter:
    ptGeom = feature.geometry()
    pt = feature.geometry().asPoint()
    nearestIds = spIndex.nearestNeighbor(pt,1)
    featureId = nearestIds[0]
    nearestIterator = lineLayer.getFeatures(QgsFeatureRequest().setFilterFid(featureId))
    nearFeature = QgsFeature()
    nearestIterator.nextFeature(nearFeature)
    shplyLineString = shapely.wkt.loads(nearFeature.geometry().exportToWkt())
    shplyPoint = shapely.wkt.loads(ptGeom.exportToWkt())
    dist = shplyLineString.distance(shplyPoint)
    shplySnapPoint = shplyLineString.interpolate(shplyLineString.project(shplyPoint))
    snapGeometry = QgsGeometry.fromWkt(shapely.wkt.dumps(shplySnapPoint))
    r = QgsRubberBand(canvas,QGis.Point)
    r.setColor(Qt.red)
    r.setToGeometry(snapGeometry,pointLayer)
    writer.addFeature(shplySnapPoint)
del writer


texts: Creating buffer based on polyline without line end caps and without dissolving the created polygons in QGIS?

texts:
I have a projected line layer (Shapefile) in QGIS. For every line segment I want to create a polygon without line end caps. The resulting polygone should not be dissolved to a single one.

I tried fTools buffering tool. This can create a buffer without dissolving, but its not possible to define the line end caps.

I tried also v.buffer of GRASS. This can create a buffer with defined line ends, but I cannot disable the dissolving.

Is there a possibility to have both: not dissolved and without caps at line end? The result should be a Shapefile with polygons (buffer) for every line segment.


texts: ArcMap toolbox failed to generate results for some users

texts:
Recently an odd problem with ArcMap happened in our office. Many people could not run the Toolbox properly. The usual map displaying was fine. But when using the Toolbox, even just a tiny job such as buffering a simple polygon, it generated no output. The Results window told that the software was running (a symbol of a hourglass under the Current Session), but it ran for a long time without any output.
Curiously, when I logged in as the administrator, it worked fine. Even when I created a new account as a standard user, the Toolbox worked. It seemed just the main user of that computer could not run the tools properly.
The data were from the local drive. We used ArcGIS 10.2.2. The system included both Windows 7 and Windows 10.


texts: Why does the map legend not update in print composer?

texts:
I'm having trouble with the QGIS map composer:

If I change the map in QGIS and reopen the map composer and click on update, the legend elements DON'T update! For example: Before one certain country was filled with blue, now with violet. The little box in the legend showing the color of that map element doesn't change. It stays blue.

How can I solve this problem?


texts: How to translate 1980 US Census Place Code to post 1990 Place FIPS?

texts:
US Census PLaces are coded differently pre-1990 and post-1990, after the introduction of FIPS Code.

For instance, for Albertville, AL:

Pre-1990 : Place Code - PLACEA  : 0030

Post-1990: Place FIPS - PLACEFP : 00988

Is there an easy way to find a correspondance between these two codes?


texts: Need unique id and to select for polygons by size after ST_Dump, ST_Intersection

texts:
I have two polygon tables, 'aa_poly_edited' and 'fields'. I am reducing the size of both (ST_Buffer), then intersecting, then using ST_Dump because after the intersection I am getting a couple of multipolygons which I don't want. Here is my query:

SELECT aa_poly_edited.id, (ST_Dump(ST_Intersection(ST_Buffer(aa_poly_edited.geom26912, -5),
  ST_Buffer(fields.geom26912, -30.48)))).geom::geometry (POLYGON, 26912) as simplegeom
FROM production.aa_poly_edited, production.fields
WHERE fields.field = 'sec35middle' AND ST_Area(ST_Buffer(aa_poly_edited.geom26912, -5)) >= 200
  AND ST_Intersects(aa_poly_edited.geom26912, fields.geom26912);

I need to add two aspects to this query and can't figure them out:

the id that is returned is no longer unique so I can't display results in QGIS. I've tried row_number () OVER () as newid but this returns ids that are the row numbers from the original feature, pre - ST_Dump so it is not unique. I also found this blog post http://longwayaround.org.uk/notes/disaggregate-multilinestrings-using-st_dump-postgis/, which uses COALESCE. There is also a SE question to this effect. Trouble is, I am new to PostGIS and SQL and I don't yet have the grasp of subqueries.
AFTER ST_Dump, I want to select only those polygons that have ST_Area >= 200  and I'm not sure where to put that in my query as well.


texts: How to add features to a PostGIS layer using Quantum GIS?

texts:
I am using Quantum GIS and the PostGIS dbms. I am able to call the PostGIS layer into QGIS but i am unable to add feature to it. I am using the method of adding features to a Shapefile but still it is not getting reflected in the layer, and I don't know how to save it to the PostGIS database if I am able to add the feature to the PostGIS layer.

I may not have used some technical term for this problem because i don't know what to call it, I am expressing the problem in my own terms.


texts: Resources for Arcpy Field Mappings String Syntax

texts:
I need to edit an existing python script to add an additional field to a field mapping function. This field mapping function does NOT use the arcpy field mapping object--it instead uses a structured string such as:

PRINMERCD \"Principal Meridian Code\" true true false 2 Text 0 0 , First,#,{0},PLSSID_1,2,3;
        PRINMER \"Principal Meridian Text\" true true false 40 Text 0 0 , First,#,{0},PrincipalMeridian_1,-1,-1;
        TWNSHPNO \"Township Number\" true true false 3 Text 0 0 , First,#,{0},TownshipNumber_1,-1,-1;
        TWNSHPFRAC \"Township Fraction\" true true false 1 Text 0 0 , First,#,{0},TownshipFraction_1,-1,-1;
        TWNSHPDIR \"Township Direction\" true true false 1 Text 0 0 , First,#,{0},TownshipDirection_1,-1,-1;

I'm able to deduce much of what I see in here, but I really need documentation. So far, I haven't been able to locate any. ESRI keeps pointing to the field mappings object, which isn't much help since it isn't feasible for me to refactor tons of code.

This article hints at my problem but then glazes over the string bit.

In summation, all I need is one document/resource explaining the syntax of this string so that I can manipulate it to suit my needs.


texts: Shapefile write protect not using Microsoft?

texts: I need to send a file to a client which can't be edited or exported. It should be write protected using a password. Is there any way to do this with a shapefile?


texts: ArcGIS Address Locator Styles & Interpolation

texts:
I'm using ArcGIS 10.3 to build an address locator service based off an address point feature class.  I want to use the Suggest operation to return possible addresses to an end user.  However, I only want ACTUAL addresses from the feature class to return, not interpolated addresses that may or may not exist.

Which ArcGIS locator styles will prevent interpolation from occurring?  Alternatively, is there a way to turn it off on any locator style?


texts: Validating feature class name in ArcGIS ModelBuilder?

texts:
I have a model which creates feature classes based on the value in an attribute (using Iterate Feature Selection). It works fine as long the field values don't have spaces, punctuation marks, etc. but fails otherwise.

I know in a script I can use arcpy.ValidateTableName() to convert "C:\Big Kalzas!" into a name that can be saved, "C__Big_Kalzas_". How does one do the same in a model?

I tried Calculate Value with an expression as simple as arcpy.ValidateTableName("%Field Name%") do defining and calling a function, but it keeps evaluating to 1.


texts: Isochrone mapping with different speeds using ArcMap Network Analyst

texts:
I'm trying to create an isochrone map using ArcMap Network Analyst (Service Area). The catch is that I want it to be travel time accurate according to each class of road's allowed speed
What would be the best way to represent my isochrone map while taking the drive time/speed into consideration?


texts: Using ArcGIS Getis Ord Gi* and IDW tools with "hotspots" where unmeasured areas are "coldspots"?

texts:
I have two datasets that I would like to analyze together. The first dataset has bird nest locations in neighborhoods (nest). The second dataset has bird feeder locations in neighborhoods (food). At each X,Y location of food I also have information about the amount of food (extent of the resource). Each of these are point shapefiles in ArcGIS.

I would like to create a smoothed surface (food raster) of food so I can extract a food index value for the location of each nest in my dataset. The food index should depend on the nearness of a birdfeeder to that nest, and the extent of food at that birdfeeder.

I used Getis Ord Gi* to generate zscores for each of my food locations (weighting each location based on extent of resource), and then used these values in IDW analysis. However, the Getis Ord Gi* zscores are both positive and negative because they compare the food value at each location to the mean food value at all of the feeders within the neighborhood. The result is that some of my feeders turn out to be "coldspots" in the resulting IDW raster.

What I'd like is a food raster where all of the food locations are "hot" (to varying degrees dependent on extent) and unmeasured locations are by definition "colder" than these (because I know there are not feeders there).

Does anyone know of a way to do this (e.g. constrain the Zscores so they are all positive or otherwise specify that the measured points are all "hot")?

Or do I need to use an alternate analysis altogether?

I have tested KDE on these data as well but would prefer IDW if possible because multiple fields of KDE are chosen at the "whim" of the researcher. I think I may have a difficult time defending choices of bandwidth, etc. with my datasets.


texts: Semi-Automatic Classification plugin: Cloud Masking errror

texts:
I'm doing Tutorial 2 out of the manual for this plugin as a GIS course lab exercise and I'm hitting an error when attempting to run the cloud masking function.

I'm using 3 bandsets of Sentinel-2 data (Bands 02,03,04,05,06,07,08,8A,11,12) of the same location on different dates. The cloud mask was converted to raster, but when running the cloud masking function I get the error below.

Error message: Warning [28]: Please define band sets with matching number of bands.

The cloud raster came from the same date bandset (bandset 1), there's no extraneous files (all .tiff format), the wavelength settings are correct for Sentinel-2, and the cloud masking page has the correct input band set selected.

QGIS version: 3.10.1
SNAP version: 6.4.2
Operating System: Windows 10 x64


texts: Masking the "infinity" part of Voronoi polygons in QGIS?

texts: I have a shapely file with several polygons for a given country, together with information on its cities. I run the Voronoi algorithm on the cities' capitals and get the polygons for it. However, polygons run off to "infinity" (outside of the external boundary of the country's polygon). How do I "mask" (is this the right word) the country so that only the Voronoi polygons are shown? In other words, how do I make polygons outside the country's boundary disappear?


texts: Running ogr2ogr script batchfile from SQL fails due to temporary drive mapping

texts:
I am trying to run the SQL procedure

EXEC xp_CMDShell 'blahblah.bat'

which needs to run an ogr2ogr script to create a kml from sql. The batch file works if I run it alone but with SQL permissions it fails. I can only think that because the first line below

pushd \\server\folder
ogr2ogr -f "KML" "\\server\folder\blahblah.kml" -s_srs "EPSG:28355" "MSSQL:server=blah;database=Blah;trusted_connection=yes" -sql "SELECT blah FROM blah" -t_srs "EPSG:28355"
popd \\server\folder

creates temporary drive mapping to the folder SQL can't find the DLL in that location.

ERROR 1: Can't load requested DLL: Y:\ogr2ogr\gdalplugins\ogr_MSSQLSpatial.dll
126: The specified module could not be found.

I've tried running the script using a SQL credential which links to a domain user which has full access to the actual folder this DLL lives in but I still get the error.

Has anyone come across this and been able to get around it? I'm assuming you can't insert the script directly into SQL.


texts: How to set background color of label using PyQGIS 3?

texts:
suppose I have a shapefile that have an attribute, and I want that attribute to be shown with white background behind the text in the viewport.

This is what I want to accomplish:

I've searched the internet and tried this: https://gis.stackexchange.com/a/277111/51131 but it only shows the labels without the background.


texts:
Is there a layout view in Quantum GIS (something similar to ArcGis layout view) ?

Anyway I am looking for tutorial or some good practices how i can prepare a nice map to print in Quantum GIS.


texts: How to remove all selection from all registered layers using QGIS plugin?

texts:
I want to remove all selection from registered layer in QGIS using Python. I tried to clear the map canvas using:

self.iface.mapCanvas().clear()
self.iface.mapCanvas().refresh()

but the selection persists.


texts: Data Driven pages, printing from mapsets not working, ArcGIS Desktop 10.2.1

texts:
Well, it's almost a year on from when I first asked this question on the ESRI user "forum" -

Data driven pages, ArcDesktop 10.2.1 - is there a limitation on how many pages can be printed?

ESRI tech support still can't find a solution.

I have almost 2500 pages to print, from several map sets built with data driven pages, all of which contain large (10x15") raster areas. So far, the "work-around" is to print 10 at a time ... and hope the printer doesn't miss a page.

For example, one mapbook MXD contains 60 data driven pages. If I attempt to print more than 10 pages, ArcMap will chew on them and appear to be processing them ... then nothing - no error message, no file sent to the printer, nothing.

With ArcMap 9.3.1, we printed MXDs containing 90 pages in one printer session (90 pages at a time), using the DSMapbook tool which was part of the ESRI Developer Samples collection. Those MXDs contained SID rasters which had been lightened, & they printed at higher quality settings.

My non-programmer's guess is that the data driven pages function is spooling multiple pages to temporary files on the computer before sending the data to the printer, and encountering some kind of limit.

After more hours with ESRI specialist printing phone tech support, we tried:

changing Windows virtual memory to 12gigs
changing the ArcGIS metafile size to 200, and back to 32
changing the default printer from a multi-function Canon to a wide-format printer-only Canon
changing driver in Print window to ArcPress
watching temp files as they were generated and discarded
watching files spool in the printer window (control panel/printers)
dropping the image resampling from 2 to 3 in the Print window

But we could only get 15 out of 60 of the 11x17" pages to print before the print process just quit with no error messages.

Have any of you good folks run into this issue? (I know I can export to PDF & print from there, but would prefer not to.)