texts: Add time dimension to vector layer in geoserver REST

texts:
I'm trying to create a layer from a PostGIS table with REST api, I can create the layer but couldn't add the time dimension.

I can configure the dimension in the interface and I did try to get the feature description but no dimensions infomation:

<layer>
    <name>focos</name>
    <type>VECTOR</type>
    <defaultStyle>
        <name>point</name>
        <atom:link xmlns:atom="http://www.w3.org/2005/Atom" rel="alternate" href="http://localhost:8080/geoserver/rest/styles/point.xml" type="application/xml"/>
    </defaultStyle>
    <resource class="featureType">
        <name>focos</name>
        <atom:link xmlns:atom="http://www.w3.org/2005/Atom" rel="alternate" href="http://localhost:8080/geoserver/rest/workspaces/terrama/datastores/view_focos13postgisstore/featuretypes/focos.xml" type="application/xml"/>
    </resource>
    <attribution>
        <logoWidth>0</logoWidth>
        <logoHeight>0</logoHeight>
    </attribution>
</layer>

Solution:
This is the xml model I'm using to do this.

<featureType>
    <name>layer_name</name>
    <metadata>
        </entry>
        <entry key="time">
            <dimensionInfo>
                <enabled>true</enabled>
                <attribute>time_attribute</attribute>
                <presentation>CONTINUOUS_INTERVAL</presentation>
                <units>ISO8601</units>
                <defaultValue>
                    <strategy>MAXIMUM</strategy>
                </defaultValue>
            </dimensionInfo>
        </entry>
    </metadata>
</featureType>


texts: Making layout size 36" x 72" in ArcGIS Desktop?

texts:
I need to make a poster to present at a conference and I'm wondering if there's anyway to do this all in ArcGIS 10.1.

Is it possible to change the layout/page size to 36x72 inches?


texts:
I have a question about GIS terminology - what is the difference between grid and raster and grid and lattice - in theory??

I thought, that raster is regular network and value is in the cell and the grid have its value in node. In this case I don´t understand what meaning have a lattice? Is the same as a grid?

Or grid is raster with value in one of the corner??

What is more confusing are some pictures, which describes difference between grid and lattice - http://www.geog.ubc.ca/courses/geob370/notes/interpolation/grid_vs_lattice.png or http://www.innovativegis.com/basis/mapanalysis/Topic18/Topic18_files/image005.jpg

In case of MATLAB example: http://www.mathworks.com/help/releases/R2013b/matlab/ref/meshgrid_ex.gif grid looks like lattice!


texts: How to implement feature search in OpenLayers with Geoserver?

texts:
I'm looking to build an application similar to this ArcGIS Server JS example which does a find task against a layer in the map. However i want to do this with OpenLayers, Geoserver and PostGIS. Some pointers would be great.

Regards

Med


texts: Why raster stack function changes the number of layers


texts: Differences between CKAN and GeoNode in spatial data infrastructures

texts:
We are planning to develop a spatial data infrastructure with some webGIS tools integrated. It would serve much vector and satellite data in numerous users. Until now, we are planning to create a project based on GeoNode because we are familiar with the platform and its components (Django, Geoserver, Postgresql, Pycsw etc). But during our research in large open data projects we realize that CKAN dominates as open data platform. Also last years it has been developed many extensions for geospatial world. So we are in a dilemma.

What are the differences between them and which is more suitable for large geospatial projects ?

At one hand we have a data management system (with spatial extensions)(CKAN) and on the other hand, a geospatial content management system (a complete package of geospatial components)(GeoNode).

I know that CKAN can harvest from GeoNode and we can use both but at first I'd like to learn about the differences between them.


texts: Can't save persisted IGeometry object into an ObjectStream

texts:
I was able to store a serialized IGeometry object into an string field in a datatable from a File Geodatabase.

For performance issues and since the object implements IPersistStream, I want to save the original IGeometry object into a blob field.

So far I'm allowed to cast the geometry to IPersistStream but haven't been able to save it to an ObjectStream

I'm using ArcObjects SDK 10.3.1. Here's the code:

// Insert Geometry to DataTable.
fieldIndex = errorsTable.FindField(m_dataBaseConstants.DataTables.Fields.Geometry);
IObjectStream objectStream = new ObjectStream();
ESRI.ArcGIS.Geometry.IGeometry geom = dataToken.Geometry;
IPersistStream persistStream = (IPersistStream)geom;
// This line throws COMException (0x80004005)
persistStream.Save(objectStream, 0);
errorRow.Value[fieldIndex] = objectStream;

I have no clue why I'm not being able to save. Can you see the problem there?


texts: Calculating new field to represent sum of values in existing field in ArcMap?

texts:
I need a way of calculating the total sum of the values in an existing field and displaying this value in the corresponding cells of another field. In other words, I want this second field to display the same value in each row, being the total of the values of the source field. I am using ModelBuilder and will be operating this function on multiple datasets, as such, this process needs to be automatic. I have no experience with coding and very little in how to implement it.

I recently asked a similar question, which was: "How to calculate the cumulative values of a field". This question was answered and resolved with the following python parser in the code block of the field calculator:

total = 0
def cumsum(inc):
 global total
 total+=inc
 return total

Where the code block was called using:

cumsum(!field_name!)

As I was able to successfully implement this parser, I would prefer if some one could adapt this example to calculate the total sum of the values in a field rather than their cumulative sum.


texts:
I am trying to upload a shapefile to ArcGIS Online, and am consistently getting an error:

"Unable to import this shapefile. This shapefile is too big to add to the map."

. I understand that there are some limitations on what can be uploaded into ArcGIS Online, namely that shapefiles must be smaller than 10MB and contain fewer than 1,000 records.

The shapefile I am uploading has 407 records. Uncompressed, all of the contingent files come to 5.59 MB. Once compressed, the zip file is only 3.55 MB (so well below the 10MB limit).

Any idea as to why I am still getting this error? Any known workarounds?


texts:
I ingested two Landsat images into PostGIS database using :

raster2pgsql -I -C -s 28356 -M -t 100x100 LS7_20150425_B40.tif public.red | psql -d postgres -U postgres
raster2pgsql -I -C -s 28356 -M -t 100x100 LS7_20150425_B50.tif public.nir | psql -d postgres -U postgres

Ingestion is going fast : around 20 seconds for each insertion.
When I did the same without tiling (-t) and vacuum (-M), it took about the same duration.

Then I want to compute a NDVI raster table from PostGIS using the following SQL command :

CREATE TABLE public.ndvi AS
SELECT ST_MapAlgebraExpr(a.rast, 1, b.rast, 1, '([rast1] - [rast2]) / ([rast1] + [rast2])::float', '32BF') AS rast
FROM public.nir  a, public.red b;

If I am not using the tiling/vacuum options at the ingestion, the NDVI raster table is computed in about 3 minutes. With the tiling and vacuum options, the NDVI raster table is computed in about one hour ! Am I missing something here ?

I decided to use the tiling option to improve the performances when displaying rasters in QGIS but if the tiling slows things down when doing some calculations there is no interest anymore to do it this way.


texts: OpenLayers override class or method in Android 2.3

texts:
I'm implementing my own image loading mechanism for OpenLayers, because lack of functions in Android 2.3 webview. It will get cached images from Android db.

I want to override OpenLayers.Tile.Image.initImage(). What options do I have?

I have try Utils for extending TMS layer, but can't extend initImage.


texts:
I'm trying to use GDAL minidriver to read tiles folders generated from Tilecache seeding.

I have set this in my GDAL source <ServerUrl>http://myserver.com/test/${z}/${x}/${y}</ServerUrl>, and when I try to use gdal_translate with this minidriver file I have 404 error:

0ERROR 1: GDALWMS: Unable to download block 0, 0.
URL: http://myserver.com/test/8/0/0
HTTP status code: 404, error: (null).

The error is right because I have another kind of structure. In example I have http://myserver.com/test/08/000/000/.

It seems to me that the output of tilecache seeding process does not have a standard TMS structure, and for this reason I have these kind of errors. Am I wrong?

My source SRS is EPSG:3004. Could it be this the reason of a non standard TMS folders structure?


texts: QGIS: separately adjust a polygon layer's fill and boundary opacity?

texts:
At 3.10.3, adjusting a polygon layer's opacity affects both the fill and the outer boundary stroke simultaneously, as shown in this screenshot:

I want to separately adjust the fill and boundary opacity.  The only way I've been able to accomplish this is to draw the layer twice, once for the fill and once for the boundary.

Is it possible to separately adjust the fill and boundary within a single layer?


texts:
I'm writing here to ask some questions for a thesis project. My ideas are confused, but I'll try to be clear.

I have a GIS project that represents some historical goods in a region. Each good has some attributes that describes the entity (like a date, a role, etc.). My task is to create a list of queries (for example: select and represents all goods that have a defined role in a defined time) and create an interface, using QGIS, that allows to represent:

A menu to choose the query,
A form to input informations,
A way to customize the presentations of results into the map and a table.

I know that my question is too broad, but to be more clear I need:

To know which tools allow to customize QGIS's interface (icons, fonts, colors, all that is possible),
Examples of QGIS projects where these two features are present: 1) customized interface, 2) possibility to query.

To help I'll attach a sketch of how the interface should be.

I'm asking if there's a way to presents the results of a query in a view like the one in my sketch, where there's the map, but also a table.


texts:
Does anyone know how to change the extents on the overview map in Flex Viewer 2.5?
I have a basemap of a single country and my overview map still wants to show me the world.
Thanks


texts: Grass GIS: ArcGIS' "Raster to polygon" equivalent

texts:
I've got a small clipped extent of a LANDSAT image and I need to have the single pixels as polygons in order to tell the overlap proportions of overlaying polygons.
In ArcGIS I used "Raster to polygon" without simplyfying polygons and used the "identity" tool then.

What would be a good equivalent in Grass GIS?

For "raster to polygon", I found "r.to.vect", but Grass always complains that the raster is "not properly thinned". I didn't understand what 'r.thin' does and I don't want ANY changes in the raster values.
Can anybody help with a "clean" solution?

EDIT:
According to Micha's answer, here's an illustration of what my final goal is:


texts: Where can I configure Geoserver getCapabilities tag ProviderSite?

texts:
I'm a beginner in GIS world and I'm having trouble with an API that reads capabilities and shows an error when tries to read the ProviderSite tag from a Geoserver getCapabilities. So, I would like to know where can I configure this tag, under ServiceProvider tag. Here is that part of the capabilities I'm serving:

<ows:ServiceProvider>
    <ows:ProviderName>http://geoserver.org</ows:ProviderName>
    <ows:ServiceContact>
        <ows:IndividualName>Whatever</ows:IndividualName>
        <ows:PositionName>Whatever</ows:PositionName>
        <ows:ContactInfo>...</ows:ContactInfo>
    </ows:ServiceContact>
</ows:ServiceProvider>

The ProviderSite should be under ProviderName. Thanks for helping.


texts: Render only features located around the the mouse cursor in QGIS?

texts:
When hovering over a vector layer in QGIS for feature selection, it would be useful, if only those in a buffer of a certain distance to the mouse cursor are displayed. Any ideas on if/how this might be available or quickly realizable?

Here is a visual example:

If the green point would be the mouse cursor, I would like to see/highlight the yellow points only.


texts: How to chessboard segment landsat pixel edges in python?

texts:
Does anyone know how to do a chessboard segmentation (in the same way that ecognition does) that will select all the edges of a Landsat pixel in python/a Jupyter notebook? Given that I can use the notebook to retrieve the imagery as a multi-band xarray, does anyone have a good recipe using scikit-image or similar to take the imagery array and run a chessboard segmentation on it, outputting the pixel edges as a shapefile? I want every pixel edge as opposed to areas of similarity.

I can do this in ecognition, and create a workflow. I'm looking for a python equivalent to the chessboard segmentation workflow in ecognition that I can apply on a large scale for multiple times. I'm not doing it in ecognition because I need to do large amounts of timesteps and large amounts of spatial area, and it's impractical in human-time-cost to do so for thousands of images.


texts: Why does every raster calcuation I try return with NAN values?

texts:
I am trying to work with DEM rasters from the USGS using the raster calculator. I have tried to change the units from meters to feet using an expression such as:

"ned10m45111h8@1" * 3.28

However this returns NAN values for every cell.  I have also tried it without the quotes around the raster name as shown in this website: http://spatialgalaxy.net/2012/01/25/using-the-qgis-raster-calculator/

Similarly, I have tried the mask code:

("ned10m45111h8@1" <= 1328.96)*"ned10m45111h8@1"

With and without quotes and get NAN returned for all values.  Is there some default setting that I need to add, or some fundamental mistake I am making?  I'm new to QGIS, but not to GIS and could easily do this with ARCGIS and spatial analyst.


texts: How to get Individual geometry elements of Fusion Table in earth engine

texts:
I have a fusion table and it content several points. I have also created NDVI layer. Now i need to print  NDVI value of individual pixels associated with point ids(taken from fusion table). But i don't know how to do that. To track the point id of fusion table i tried to print point features. I have written the following codes for the same. Console pictures is attached to which ids i want to get.Also sharing the GEE linkCode link [Exmp: [Need to know the ndvi value for id 3]]

    // ---------Cloud Masking function

function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = ee.Number(2).pow(10).int();
  var cirrusBitMask = ee.Number(2).pow(11).int();

  // Both flags should be set to zero, indicating clear conditions
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));

  // Return the masked and scaled data.
  return image.updateMask(mask).divide(10000);
}

var boundary = ee.FeatureCollection('ft:1T-XrxMki2k49fhQdhMxesWKBoyyUJiKN70pErh9z');// Boundary
var gtpoints = ee.FeatureCollection('ft:1ZvFeFO5jHGo3VExlKNc-kdL-3hmJSeyfLzqLvHhy');//Co-ordinates
print(gtpoints);
var test = gtpoints.features;
print(test);
Map.centerObject(boundary);
// ---------  Import of Images (Sentinel _ 2 multispectral)

var image = ee.ImageCollection(sent2img
.filterDate("2018-09-22","2018-10-23")
.filterBounds(boundary)
.map(maskS2clouds)
.sort("CLOUD_COVERAGE_ASSESSMENT")
.median()
);

//----------  Preprocessing

var mosaic = image.mosaic()
var clip = mosaic.clip(boundary);
print('clip', clip);

//----------   FCC creation and visualisation of AOI
Map.addLayer(clip, {bands: ['B8','B4','B3'], min: 0, max: 0.3},'clip');

var ndvi = clip.expression(
    ' ((NIR - RED) / (NIR + RED))', {
      'NIR': clip.select('B8'),
      'RED': clip.select('B4'),
}).rename('nd');

var palette = ['196fda', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];

Map.addLayer(ndvi, {min:0, max:1, palette: palette},"NDVI");


texts: Symbolizing changes in housing condition over time using ArcMap?

texts:
I have two polygon layers, one from 2015 and the other from 2018. Each layer has an attribute field with a number from 1-14 which describes the physical conditions of houses. There are 15,000+ houses in each layer.

What I'd like to do using ArcMap is to create a map that answers the following question:

Are the housing conditions in 2018 better, the same, or worse than 2015?

Image below: YEAR 2015

Image below: YEAR 2018


texts: How to filter OSM data spatially (more complex than by polygon/bounding box)?

texts:
Tools like osmosis, osmfilter, or osmconvert allow just clipping data to a defined bounding box or a polygon (besides filtering by attributes).

What is the best way to delete all objects inside multiple defined areas? In other words: I want do do an erase operation on OSM data.
Is there a simpler way than importing OSM data to a database, processing there and writing back to a file?


texts:
I need to analyze data uploaded from vehicle GPS devices to detect travel patterns and then attempt to learn about the normal travel patterns for that vehicle so that predictions can be made about maintenance, refueling etc. Some vehicles are domestic (exhibiting normal home->work->home travel patterns) and some are commercial with scheduled routing. I have loaded the GPS lats and longs to a database, and then constructed an ST_LineString with the data points for each vehicle (approx 1k points per vehicle per day). Vehicles will never have identical travel patterns (GPS accuracy, irregular ping times etc) and the patterns will be visible as subsets of pings - for example a one hour commute between point A (home) and point B (work) between 8am and 9.30am each morning, mon-fri.

Now the difficult part - I am considering the options. I have to write this in a database procedural language using the usual toolbox of procedures, cursors, iteration etc. I could use Python if it is helpful - I'm a database guy trying to think laterally and spatially.

A travel pattern could be identified by

(i) a similar sequence of points that can differ by a tolerance (say +/- 300 meters)

(ii) parking locations - if a car is parked at the same location for 8 hours a day

(iii) temporal- the vehicle starts at trip at 0830 every day, +/- 10 minutes.

I am planning to write a procedure that checks for each of these patterns for each point (or set of n points) by comparing it with all other points (or set of points) by iterating through the pings for each day in sequence. The easiest use case is to find recurrences of a vehicle parked at the same place each day. Detecting similar movement patterns is more complex. Can anyone suggest a technique for this?

Thanks!


texts: DEM provided contains high elevations (trees) that I wish to remove - is this possible

texts:
I wish to remove trees off my 3D aerial (high elevations). The files provided are  .tif - unfortunately the trees are causing quite a bit of noise. Is there a way to set an elevation to remove them? Area is not very sloped so the trees would be chopped.

I also have a regional .tif which does not include the trees - could i stack both .tif together to somehow remove the trees...

Sorry i am fairly new to QGIS any help would be wonderful.


texts: Plot points from a Sqlite (non spatial) table as a query on the canvas instead of making a new table

texts:
I have a table in an SQLite database with a non spatial table that contains 6 sets of coordinates, representing different features within a record. I would like to generate a series of layers (lets say 6 for argument's sake) and because my dataset is going to grow I don't want to create new tables each time...

I think I'd like to create query that has a geometry column plus the related columns taken from my record. That would (I hope) allow me to requery the table when the dataset is updated...

I thought about:

Select 'Point (' || x1 || ' ' || y1 ||')' as geometry, name, description, uuid from mytable

That gives me the columns, with a WKT column masquerading as my geometry. Where I get stuck is throwing this on to the canvas...
I tried 'load as layer' in DB manager, but when I select geometry as my geometry column nothing is loaded to the canvas.


texts:
Does WMSGetFeatureInfo need a proxyhost to allow remote server accessing the data? I use ProxyHost to allow WFS requests, but seems its not working with WMSGetFeatureInfo. the local machines retrieves data from GeoServer successfully but external machines doesn't!

Here is the select feature code:

    control = new OpenLayers.Control.WMSGetFeatureInfo({
        layers: [layer],
        autoActivate: true,
        maxFeatures: 3,
        hover: false,
        infoFormat: "application/vnd.ogc.gml",
        eventListeners: {
            "getfeatureinfo": function (e) {

                var activeLayer = getActiveLayer();
                if (activeLayer != "") {
                    var param = {
                        REQUEST: "GetFeatureInfo",
                        EXCEPTIONS: " application/vnd.ogc.se_xml",
                        BBOX: map.getExtent().toBBOX(),
                        SERVICE: "WMS",
                        VERSION: "1.1.1",
                        X: e.xy.x,
                        Y: e.xy.y,
                        INFO_FORMAT: 'application/vnd.ogc.gml/3.1.1',
                        QUERY_LAYERS: activeLayer.params.LAYERS,
                        Layers: activeLayer.params.LAYERS,
                        FEATURE_COUNT: 6,
                        SRS: 'EPSG:3857',
                        STYLES: activeLayer.params.STYLES,
                        WIDTH: map.size.w,
                        HEIGHT: map.size.h
                    };

                    param.cql_filter = activeLayer.params.CQL_FILTER;

               //     if (activeLayer.params.LAYERS == 'AIS:favoriteVessels');
                //    param.Viewparams = 'UID:' + userid;

                    OpenLayers.loadURL(GeoServerWMS, param, this, setHTML, setHTML);
                    OpenLayers.Event.stop(e);
                }
            }
        }
    });

map.addControl(control);
control.activate();

and the GeoServerWMS is specified by:

<php
if ($_SERVER['SERVER_NAME'] == 'MYPUBLICIP')
  GeoServerWMS = 'http://MYPUBLICIP:8080/geoserver/wms'
else
  GeoServerWMS = 'http://privateIP:8080/geoserver/wms'

here is my allowedHosts

allowedHosts = ['localhost', 'localhost:8080','localhost:80', 'localhost:8080/geoserver',
                'PrivateIP','PrivateIP:8080', 'PrivateIP:8080/geoserver',
                'PublicIP','PublicIP:8080', 'PublicIP:8080/geoserver',
                'www.openlayers.org', 'openlayers.org',
              ]

EDIT
Even this event register is not working:

selectControl = new OpenLayers.Control.GetFeature({
    id: 'selectControl',
    protocol: OpenLayers.Protocol.WFS.fromWMSLayer(test, {
    url: GeoServerWMS,
        featureType: "layer"
    }),
    hover: false,
    click: true
});

Here is the error I get: Some unexpected error occurred. Error text was

  urlopen error [Errno 110] Connection timed out


texts: Calculating total area of three partially overlapping polygons (without counting area of overlap more than once) using ArcMap

texts:
I want to calculate the area of three 1 km buffers without counting the area of overlap more than once. I have used the merge tool to combined the three polygons and then I have used the dissolve tool.
Is this the correct way to calculate the total area of the three polygons without counting the area of overlap more than once?
I wasn't sure if there was another tool which would be better suited to this.


texts: Use multiple attributes from table to select and merge multiple shape files

texts:
I'm totally struggling with a model build in ArcGIS 10.3.

Here’s what I’d like to do.

Take a table (analysis_table). Iterate over the rows in the table by a unique grouping identifier. From each row, use the values of two different fields, which are geoids, and select shapefiles from a separate database by those geoids. Create a new merged shapefile based on the two different shapefiles from the step immediately prior.

Here’s what happens:

I am able to parse out the value of two different variables within the row. This appears in the model as get field value (2) and get field value (3). These values are then fed into “select layer by attribute”, using a %value% expression.

Unfortunately though, select layer by attribute either has me make a new selection with each pass, or if I add to selection, it grows the selection continuously. In other words, rather than select a shapefile for each row in the iteration, it does one for each value that is passed from “I_analysis_table_OBJECTID”.

I want to merge the two shapefile outputs of “county_500k” (selected from my row iteration), save them out, then clear the selection for the next row pass, but I don’t know how or if that is possible.

You are right that I wasn't clear enough. I actually used different inline variable expressions for each Get Field Value, but I was getting sequentially more shapefiles added on with each row pass in the iterator. What I wanted was a new combined shapefile for each row pass of the iterator.

Today I took a different tact, and created sub-models to link each different attribute to each county, then bringing them back into the main model. I used a combination of "new selection" and "add to selection" in the select by attribute drop down options, which allowed each pass to start with a new selection, then add the sequential counties. You'll see that I added a third sub-model, but the syntax is essentially the same as the original problem I was trying to solve with two.

If it would be helpful to others, I can give more detail on the inline variable substitutions.

I'm including screenshots of the base model and one of the sub-models.


texts:
Treat me as a beginner.

I have this method, where the bandArray is the original image double[] array

private double[] ApplyAverageSmoothingNxN(double[] bandArray, int width, int height, int n)
{
    //logic inside here (not relevant)
}

This method returns me a double[] array which is actually the raster file. After the methods returns I try to save the output using Gdal like so:

composedDs = drv.Create(path + "smoothed.tif", XSize, YSize, 1, DataType.GDT_Int16, null);
composedDs.SetGeoTransform(transform);
composedDs.WriteRaster(0, 0, band.XSize, band.YSize, result, band.XSize, band.YSize, 1, null, 0, 0, 0);
composedDs.Dispose();

Now, after I open the output tif image, I see everything black (all 0's). But if I try to run my code again, the raster has all it's values correct. However, if I try to assign all values of the output double[] array like so:

for (int i = 0; i < array.Length; i++) {
     array[i] = 10 //some value.. could be anything else
}

After I save it with Gdal.Create() Arcmap displays pixel values correctly. What is happening ?

Also my input image is a ".jp2" file while the output is done with "GTiff" driver as ".tif".


texts: "TypeError: this.options is null" in OpenLayers.Protocol.HTTP

texts:
I'm using OL 2.11 with GeoExt 1.1. After testing and modifying this GeoExt example: http://geoext.org/lib/GeoExt/widgets/form/FormPanel.html, I had problems with destroying and reading the new protocol, I'm getting "TypeError: this.options is null      HTTP.js (line 180)" in firebug. What I understand from this line 180 is that the new options (in the new protocol) are not there, not sure if this is the problem and how to solve it. This is the original protocol (created outside the FormPanel):

protocol = new OpenLayers.Protocol.HTTP({
    url: '/fs/',
    format: new OpenLayers.Format.GeoJSON({
        ignoreExtraDims: true,
        'internalProjection': new OpenLayers.Projection('EPSG:900913'),
        'externalProjection': new OpenLayers.Projection('EPSG:4326')
    })
});

This is the code at the end of the FormPanel where the new protocol is created:

    buttons: [{
        text: 'Search',
        handler: function() {

            comboLayer = Ext.getCmp('idcombo').getValue();
            keyword = Ext.getCmp('idtextfield').getRawValue();

            protocol.destroy();

            protocol = new OpenLayers.Protocol.HTTP({
                url: '/fs/' + comboLayer + '?format=GeoJSON&comments__ilike=' + keyword + '&queryable=comments',
                format: new OpenLayers.Format.GeoJSON({
                    ignoreExtraDims: true,
                    'internalProjection': new OpenLayers.Projection('EPSG:900913'),
                    'externalProjection': new OpenLayers.Projection('EPSG:4326')
                })
            });

            protocol.read();

            formPanel.search();
        }
    }]

In the case I remove "protocol.destroy();", both the original and the new protocols are sent, but always the original is sent after the new one, this happens because the request is AJAX? is it possible to send the new after the original?

I'd appreciate some support on this, thanks in advance,


texts:
I would like to use OpenLayers 3 spyglass example inside my project.

The problem is when I want to use selection for layer, which will be inside spyglass circle. I can't get it working, I don't know how to update it on every user's selection correctly.

Anyone more skilled can help me?

Here is my fiddle with it partially working:
http://jsfiddle.net/dkcz/8a3pdavv/


texts: Is it possible to have a locked/different first page with QGIS Atlas printing?

texts:
It is possible set the first page of a QGIS Atlas different to the other pages? For example I want to create an Atlas document where the first page has a text introduction and a general map (not controlled by Atlas), then from the second page and onwards I want detailed maps but in this case controlled by Atlas.

I have tried to do that setting two pages at the "Composition" tab where in the first page is the introduction and in the second page the map controlled by Atlas, but when I create the pdf document Atlas repeats the first page (page #1) every one page. It is possible to "lock" or "fix" the first page to prevent it from repeating in the rest of the document?


texts: I am new to geoNetwork, I want to know why geoNetwork use H2 database, Is there any special reason not to use mysql or pgsql? what are the special features in H2 database in GIS, I search about this in google but i didn't get clear answer.


texts: Processing LWGEOM Provider crashes QGIS 2.18

texts:
I am trying to use the LWGEOM processing algorithms in QGIS 2.18 running on OS X (El Capitan). I have installed the liblwgeom library (v. 2.1.5) using a Homebrew formula and successfully set its provider path under Processing > Options. When I run the algorithm, QGIS crashes immediately.

What might be causing this crash, and how can I confirm that I've installed the algorithm correctly?


texts: Geoprocessing tool default path access every one?

texts: I am using Geoprocessing tools i have face one problem with default gdb access in my system (default gdb means C:\Documents and Settings\\My Documents\ArcGIS\Default.gdb ) how can access different users same default gdb file. Example my user is schalla it is working fine same tool use another person it is getting error message like attachment


texts: How to find frequency of a field based on another attribute?

texts:
I'm working in ArcMap. I have a list of people and the states that they live in. Small Example:

I want to know what is the most common name in each state.

I thought I would use the frequency tool to return unique values of names, but I want to know the frequency for each state. Is there a way to do this without doing the frequency tool 50 times for each state?


texts: How can I measure area from geographic coordinates?

texts: If I have polygons in geographic coordinates (WGS84), how do I measure the total area each takes on the surface of the earth, taking into account the curvature of the earth?


texts: Error sld generated for mbstyle in GeoServer

texts:
I have activated the mbstyle plugin for GeoServer. This is part of a mbstyle, it contains a filter:

"filter": ["all", [">=", "scalerank", ["zoom"]]],

Which means I want to render only those features whose scalerank is >= the current zoom level.

And I found a sld file named the same as the mbstyle in the $geoserver_data_dir/styles/xxx.sld, this should be generated by geoserver, and the above filer is transformed as this:

  <ogc:Filter>
    <ogc:And>
      <ogc:PropertyIsGreaterThanOrEqualTo>
        <ogc:PropertyName>scalerank</ogc:PropertyName>
        <ogc:Literal>zoomLevel([env([wms_scale_denominator])], [EPSG:3857])</ogc:Literal>
      </ogc:PropertyIsGreaterThanOrEqualTo>
    </ogc:And>
  </ogc:Filter>

When I preview the layer, I got the following error:

  java.sql.SQLException: org.postgresql.util.PSQLException: ERROR: invalid input syntax for integer: "zoomLevel([env([wms_scale_denominator])], [EPSG:3857])"

How can I fix this?

Update for answer by @Andrea Aime:

I found the class ZoomLevelFunction extends FunctionExpressionImpl which provide the zoomLevel function, so the transformation in the sld is right/expected. But Geoserver pass this expression to database for features querying directly instead of  evaluating it using the function. So is it the problem of geoserver or mbstyle?

Also follow this PR mentioned by @Andrea Aime, I change my filter to :

"filter": ["all", [">=", ["get","scalerank"], ["zoom"]]],

And the sld changed to:

  <ogc:Filter>
    <ogc:And>
      <ogc:PropertyIsGreaterThanOrEqualTo>
        <ogc:PropertyName>property([scalerank])</ogc:PropertyName>
        <ogc:Literal>zoomLevel([env([wms_scale_denominator])], [EPSG:3857])</ogc:Literal>
      </ogc:PropertyIsGreaterThanOrEqualTo>
    </ogc:And>
  </ogc:Filter>

But I got error:

  This style specifies an attribute name property([scalerank]),not found in the xxx layer`.


texts: Bubbles getting shifted for WMS getMap if POST is used instead of GET, using gutter in both cases

texts:
Bubble positions are wrong in case of WMS POST request. For example, a bubble is getting shifted from 31.82, -96.84(lat-long) to 31.77, -96.77.
I am using openlayers 2.13.1.

The layer creation Code is:

   wms: {
      options: {
        gutter: 200,
        tileOptions: {maxGetUrlLength: 2048},
        ratio: 1,
        singleTile: true
      },..

I am using projection:"EPSG:3857". SingleTile and ratio options are not being used originally, but for experimental purpose, but irrespective of that results are same.

As per observation,

  GET and POST requests(using cURL) for geo-server give same images.
  If I don't use gutter, it gives correct result for POST (and GET also)
  POST url contains BBOX info, I am not sure use of that but it is different that FORM > BBOX, in case of gutter. If no gutter is used then both are equal.
  I just went through OpenLayers code. It seems while determining GET URL (which gets converted into POST FORM) gutter is used to adjust BBOX, it is not happening in IFrame, I guess.


texts: Snap a point to a line and get the line ID (using the ArcGIS REST API)

texts:
I have an external work order management system in which users create new sidewalk work orders (points).

When a user creates a sidewalk work order, a Jython script will invoke a URL that passes the X & Y coordinates of the work order to a map service/ArcGIS REST API (10.7.1).

In the map service/ArcGIS REST API, I want to process the X & Y coordinates as follows:

Find the nearest sidewalk to the X & Y coordinates (within a 5-metre distance) and get its ID.

Sidewalks are contained in a polyline EGDB feature class/feature service.

Snap the X & Y coordinates to the closest part of the sidewalk.
Return the sidewalk ID and the snapped coordinates to the work order management system.

Is it possible to do this using the ArcGIS REST API?

Unfortunately, I only have the Standard ArcGIS Server license. So, while I can use a geoprocessing service, I cannot use the Near tool.


texts: getting the correct location data with user input

texts:
I have some code others have provided, which I modified in order to run a Python script in a toolbox in ArcMap which accepts user input for a starting point, distance, and angle to create a shapefile for a line. The script works OK except the starting point (origin_x, origin_y), which I am entering in decimal degrees, is always wrong when I place the shapefile in a map dataframe. How can I correct this code so that the math can work correctly to produce the right starting point location? It seems converting origin_x and origin_y into radians is the wrong way to go about it, but I'm not sure...

import arcpy
from arcpy import env
from math import radians, sin, cos, pi

#Set workspace
env.workspace = "G:\\Delapaix\\Geog_586_FinalProject\\line_connect"

origin_x   = arcpy.GetParameterAsText(0)
origin_y  = arcpy.GetParameterAsText(1)
distance  = arcpy.GetParameterAsText(2)
angle  = arcpy.GetParameterAsText(3)

# convert latitude and longitude to radians
origin_x = float(origin_x)/(180/pi)
origin_y = float(origin_y)/(180/pi)

# calculate offsets with light trig
(disp_x, disp_y) = (float(distance) * sin(radians(float(angle))), float(distance) *  cos(radians(float(angle))))
(end_x, end_y) = (float(origin_x) + disp_x, float(origin_y) + disp_y)

output = "offsetLine.shp"
arcpy.CreateFeatureclass_management("G:\\Delapaix\\Geog_586_FinalProject\\line_connect",output, "Polyline")
cur = arcpy.InsertCursor(output)
lineArray = arcpy.Array()

# start point
start = arcpy.Point()
(start.ID, start.X, start.Y) = (1, origin_x, origin_y)
lineArray.add(start)

# end point
end = arcpy.Point()
(end.ID, end.X, end.Y) = (2, end_x, end_y)
lineArray.add(end)

# write feature to the shapefile
feat = cur.newRow()
feat.shape = lineArray
cur.insertRow(feat)

# code for buffering the line and including features in buffer zone would go here

# yes, this shouldn't really be necessary
lineArray.removeAll()


texts: Setting Data Driven Pages rotation when index has extra page polygons edited into it?

texts:
I have created an index to use with data drive pages. I used strip map index and then altered the geometry of the index polygons to better fit my project area. I also added additional index polygons. As a result I copied the geometry of the polygons and created a new index shapefile. In this new shapefile I populated the attribute table with labeling along with before and after label (since some of the maps are sequential for match lines).

I now need to add a field which calculates the data frame rotation for each polygon. Strip Map index does this automatically, but under the circumstances, I need to do this semi manually. How do I populate an attribute field with the data frame rotation angle?


texts: How to overlay and combine a choropleth and a shaded raster in Python?

texts:
I want to make a choropleth map of regions but shade it according to local population density. This is not the problem of masking a raster with a shapefile but to overlay two representations of data over the same area. One (population density) is as good as continuous for our purpose and could be rasterized in whichever object is necessary. The other data is coarse, at the level of subregions, and thus comes vectorized.

Colors (hues) could come from the choropleth (thus stay constant within subregions) but the transparency of each pixel should vary inversely with local population density, as approximated by a local average in a rastering call.

See the two independent plots below. Rasterio or the HoloViz family (GeoViz and Datashader) look relevant but I could not make them accomplish this.

The raster comes from Facebook's high-resolution maps of residence, here for Hungarian youth.

import dask.dataframe as dd
df = dd.read_csv("HUN_youth_15_24_2020-02-01.csv")
import datashader as ds
from datashader import transfer_functions as tf
from datashader.colors import Greys9
Greys9_r = list(reversed(Greys9))[:-2]
plot_width  = int(750)
plot_height = int(plot_width//1.2)
x_range, y_range =(16.141528,22.883194), (45.755139,48.582917)
cvs = ds.Canvas(plot_width=plot_width, plot_height=plot_height, x_range=x_range, y_range=y_range)
agg = cvs.points(df, 'longitude', 'latitude',  ds.mean('population'))
img = tf.shade(agg,cmap=Greys9_r, how='log')

And this produces a choropleth map of the relevant subregions of Hungary (here with mock data), using the shapefile from GADM 3.6:

import pandas as pd
import geopandas as gpd
import xarray as xr
import numpy as np
import mkl_random
shp = gpd.read_file('gadm36_HUN_2.shp')
GID_2 = shp['GID_2']
values = mkl_random.uniform(0.5, 20, 168)
arr = xr.DataArray(values, dims='GID_2', coords={'GID_2':GID_2})
shp = shp.set_index('GID_2')
shp = shp.assign(myvar=arr.to_series())
shp.plot('myvar')


texts: NIR, R, G, B multiband (4 Band) GeoTIFF gdal_translate colorinterp issues

texts:
I originally create a 5 band image (NIR, R, G, B, NDVI) from Sentinel 2 Sat imagery that has valid data in the bands to work with. I'm stacking them together with rio stack command.

rio info 5band.tif

{"count": 5, "crs": "EPSG:32614", "colorinterp": ["gray", "undefined",
"undefined", "undefined", "undefined"], "interleave": "pixel", "dtype":
"uint8", "driver": "GTiff", etc...}

NIR, G, B in QGIS from 5 band:

However, if I try to use either rio cli or for instance the following gdal_translate command to create a new GeoTIFF with only the first 4 bands (NIR, R, G, B), the colorinterp changes for the bands and the image output appears to be substantially lighter because its setting the 4th band as alpha on the colorinterp which gets loaded as transparency.

gdal_translate -b 1 -b 2 -b 3 -b 4 5band.tif 4band.tif

{"count": 4, "crs": "EPSG:32614", "colorinterp": ["gray", "green",
"blue", "alpha"], "interleave": "pixel", "dtype": "uint8", "driver":
"GTiff", etc...]}

NIR, G, B in QGIS from 4 band:

Any ideas why this is happening and what I can do to prevent it from registering the 4th band as alpha?


texts: ArcGIS Pro Arcade expression to Maintain line-to-point connectivity

texts:
A customer was using a geometric network to maintain connectivity between lines and points. Now they are using ArcGIS Pro, which does not support geometric networks.

They cannot use the Utility Network since they don't have an ArcGIS Enterprise license.

We are trying to handle connectivity issues with Arcade expressions in attribute rules.

Using this approach, when we move a point feature, is there a way to make the point that's connected to it move also - the same way it does in a geometric network?


texts: How to project Lambert Conformal Conic to WGS84 in Python?

texts:
I am using pyproj to re-project Canada's Census Tract boundary shapefiles that are currently in "PCS Lambert Conformal Conic" system. I want to convert them to latitude and longitude WGS84 values.

Is "PCS Lambert Conformal Conic" the same as "Canada Lambert Conformal Conic" system (since I could not find a code for the former)?
If yes, the code I found for "Canada Lambert Conformal Conic" is ESRI:102002 and for WGS84, it is EPSG:4326.

My coordinates in Lambert system are x1, y1 = (8411012.879999999, 1438599.017145). After reprojection using the above codes, I get (-24.41452149523307, 11.570599554812532) which is definitely not in Canada. I get the same results when I use these codes in pyproj as well as on this website: https://epsg.io/. Since the coordinates are from Halifax, I expected the reprojected coorindates to be something near 44.6488° N, 63.5752° W (Halifax).
What am I doing wrong here?


texts:
I can toggle the active basemap programmatically using the code below, but how can I also include a widget to show or hide the basemap labels programmatically?

myApp.widgets.basemapManager.basemapGallery.activeBasemap = {
  portalItem: {
      id: baseMapId
    }

I see no labels object or properties in the dev console for activeBasemap, how can I do this?

I am using ArcGIS JS API 4


texts: Getting spatially referenced tiles from WMTS

texts:
I have been able to download tiles from a WMTS request in .tif format using owslib in standalone python:

from owslib.wmts import WebMapTileService

wmts = WebMapTileService('wmts_url')
img = wmts.gettile(layer='v0Lw1953KMdH',
                   tilematrixset='WGS84',
                   tilematrix='16',
                   row='20480',
                   column='32001',
                   format='image/tif')

out = open('C:\\out\\test_out.tif', 'wb')
bytes_written = out.write(img.read())
out.close()

This .tif file is not spatially referenced. The WMTS Layer has the tag:

<ows:WGS84BoundingBox crs="urn:ogc:def:crs:OGC:2:84">
    <ows:LowerCorner>-7.793077 49.852539</ows:LowerCorner>
    <ows:UpperCorner>1.790425 60.894042</ows:UpperCorner>
</ows:WGS84BoundingBox>

and the WMTS TileMatrixSet has the tag:

<ows:Title>WGS84</ows:Title>
<ows:Abstract>WGS84 EPSG:4326</ows:Abstract>
<ows:Identifier>WGS84</ows:Identifier>
<ows:SupportedCRS>urn:ogc:def:crs:EPSG::4326</ows:SupportedCRS>

which to someone of my limited knoweledge indicates that there could be a way to obtain a spatially referenced .tif from these tiles. However, I haven't been able to obtain a spatially referenced .tif, despite trying to include the crs in gettile(). Where am I going wrong?

Failing this, are there any other methods that I could use to obtain spatially referenced tiles?


texts:
I'm trying to render a KML layer (with timezones) in Google Maps API v3. The KML layer works fine in 3D (in Google Earth), but when I render it in 2D, the timezone lines that meet at the North Pole are spread out in a really weird way:

Here is the correct rendering in 3D:

How do I get fix that weird spread?


texts: Can shapefiles be converted to .TAB files?

texts:
I read this response to this question, "Yes, if you have access to MapInfo, you can directly read Shapefiles (as .shp is an "open" format) and save it as TAB files!"

Once it is saved as a .TAB, are you then free to modify the data in MapInfo?


texts: How to find out whether a lat/long position resides within a Polygon?

texts:
I'm working on a system that stores latitude/longitude information of addresses in a PostGIS table. To store Lat/Long in geometry data type I'm using ST_GeometryFromText function.

For example the following function call gets the geometric value against the specified Long/Lat position:

 myPointGeo =  ST_GeometryFromText('POINT(40.758871 -73.985114)');

Similarly, I convert a Polygon into geometric representation as follows:

myPolygonGeo =  ST_GeometryFromText('POLYGON ((40.7566484549725 -73.9878561496734, 40.7556894646734 -73.9853026866913, 40.7545841705587 -73.9860537052154, 40.7548036054111 -73.9881458282471, 40.7559820394514 -73.9887895584106, 40.7566484549725 -73.9878561496734 ))')

I want to find out whether the above lat/long position resides within this Polygon or not. For this I'm using ST_Within function. But though 'myPointGeo' actually resides within the Polygon, ST_Within is returning false. I'm using ST_Within this way:

 St_Within(myPointGeo,myPolygonGeo)

What am I doing wrong here? Should I use some other function for this purpose?


texts: Loading shapefiles for reprojection using sf

texts:
I'm trying to automate my workflow, but having trouble with one step. My data come as processed shapefiles saved in .RDA files; I load them to a PostGIS system, cleaning up little issues (spTransform() to a local CRS, as.numeric() on a key field, etc).

When I simply load(c:/data/regionX-yearY.rda) they load as SpatialPolygonsDataFrame so I can simply
regionX-yearY_3035 <- spTransform(regionX-yearY, "+init=epsg:3035")
... but for some reason looping over files from list.files() causes them to be loaded as characters, or their filehandle or something. The sf package errs with

  "unable to find an inherited method for function ‘spTransform’ for
  signature ‘"character", "character"’"

Ideas as to what's going on here? Obviously I took some sort of implicit shortcut by hand.


texts: Very slow loading while fetching WMS layer as image/svg ol5

texts:
I want to fetch WMS layer as SVG image but it is loading very slow. How to resolve this issue?

Here is my code

var grid_10hact = new ol.layer.Image({
  title: '10 Hactare Grid',

    source: new ol.source.ImageWMS({

        url: 'http://localhost:8080/geoserver/wms',

        params: {
            'LAYERS': 'soil_health:CG_STATE_GRIDE_10HA_WGS84',

            'VERSION': '1.1.1',
            'FORMAT': 'image/svg',

        }
    }),
    opacity: 0.6
});


texts: Standard name for a latitude/longitude pair

texts:
Is there a standard name for a latitude/longitude pair? I ask because a lot of mapping APIs (e.g., Google Maps, Mapbox, MapQuest) represent them using data structures named LatLng. That naming is very to-the-point but it just seems lazy.

Meanwhile, other libraries (e.g., OSM, Proj.4, JTS, ArcGIS, PostGIS, Apple Maps) refer to them as "coordinates" or "points." Is there one particular term that is more accurate than the rest for describing a point on the globe that is created from/represented by a latitude and longitude?


texts:
I am trying to understand how the extent works in the javascript api.  I create an extent like the following:

var beginExt = new esri.geometry.Extent(-159.6716766, 19.520, -67.843, 64.8348, new esri.SpatialReference({wkid:4326}));
var revExt = esri.geometry.geographicToWebMercator(beginExt);

I add it to the map but when I display the current extent with this:

function showExtent(extent) {
        var s = "";
        s = "XMin: "+ extent.xmin.toFixed(2) + " "
           +"YMin: " + extent.ymin.toFixed(2) + " "
           +"XMax: " + extent.xmax.toFixed(2) + " "
           +"YMax: " + extent.ymax.toFixed(2);
        dojo.byId("info").innerHTML = s;
      }

I get this as output:

  XMin: -9384566.21 YMin: 4732448.76 XMax: -7756763.25 YMax: 5357397.90

This is not even close.  It appears to be a different format maybe using a different spatial reference?  I would like to to be able to see the extent so I can set the extent I need for the map.  However when I replace my existing Extent with these output numbers the page does not display.  How can I get my extent in the correct format so I can see an extent and use those numbers to set it for my map?


texts: Using a join in Arc with feature classes containing a definition query always return an inner join

texts: At Arc 10.4, I am attempting to do an outer join (using the "Keep all records" option) between two file geodatabase feature class layer tables.  Each of the two tables has an identical Definition Query applied to it. The join completes without an error statement from Arc.  However, the resulting join is an inner join (containing only the matching records), instead of the expected outer join. I conducted an internet search, which did not return anything useful.  Is this a known Arc bug?


texts:
I have 2 layers ( geometric type - point), let's name it layer A and layer B. I need to find which points from layer B don't have their equivalents points from layer A. So i`ve made 10 meter buffer around points A and found which interestcs with points from layer B. So now I have a group of points from layer B which are intersected with buffers so they have equivalents points in A, and I can also find points from layer B which don't have equivalent points in A. But there is a chance that 2 buffers from layer A can intersect with 1 (the same) point from layer B. How can I find which 2 buffers from layer A are intersected to the same point. Or maybe on the other way - which point from layer B is intersected to 2 or more buffers?

The problem is in the 3 example, where 2 buffers are intersected to one point. How i can find that buffers?

I am using QGIS.


texts: QGIS difference output has same number of features

texts:
Having a problem with the QGIS Difference feature. Two inputs: a point and a polygon file. I am trying to generate a list of the points that are not covered by the polygons.
I am using the "use only selected features" option on the point shapefile, as I want only points that are within the US. When I use "only selected", the difference output looks correct on the map but my attribute table contains all selected entries in the point file!

When I do not use the "selected features" option, the attribute table contains only the features that are outside of the polygons (it works correctly). I can then export and filter the list appropriately.

Question:
What am I doing wrong? (if anything)? Why does the map show the correct result, but with entries in the attribute table that are not mapped? Is this a bug that I should report? (Using 2.12.3-Lyon on Win7)


texts:
I am wondering where i can  find the map service address for basemap "street" found at \viewer\js\config\viewer.js ? (Configurable Map Service v 1.3.4)

>     define([    'esri/units',    'esri/geometry/Extent',    'esri/config', 'esri/tasks/GeometryService',
> 'esri/layers/ImageParameters' ], function (units, Extent, esriConfig,
> GeometryService, ImageParameters) {
>
>         // url to your proxy page, must be on same machine hosting you app. See proxy folder for readme.
>         esriConfig.defaults.io.proxyUrl = 'proxy/proxy.ashx';
>         esriConfig.defaults.io.alwaysUseProxy = false;
>         // url to your geometry server.
>         esriConfig.defaults.geometryService = new GeometryService('http://tasks.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer');
>
>         //image parameters for dynamic services, set to png32 for higher quality exports.
>         var imageParameters = new ImageParameters();
>         imageParameters.format = 'png32';
>
>         return {
>             // used for debugging your app
>             isDebug: true,
>
>             //default mapClick mode, mapClickMode lets widgets know what mode the map is in to avoid multipult map click actions from
> taking place (ie identify while drawing).
>             defaultMapClickMode: 'identify',
>             // map options, passed to map constructor. see: https://developers.arcgis.com/javascript/jsapi/map-amd.html#map1
>             mapOptions: {

        **basemap: 'streets',**  //where is it defined ?

>                 center: [-99.59179687497497, 139.09596293629694],
>                 zoom: 5,
>                 sliderStyle: 'small'
>             },


texts: Performing join on simultaneous spatial and attribute criteria in ArcGIS for Desktop in Modelbuilder?

texts: A similar question was posted here, but I have an additional problem that seems to make the answers there unlikely to work.  I am interested in performing a spatial join between point feature class A (target) and point feature class B in order to get the feature in B that is closest to each feature in A.  However, I would like it to be the case that a feature in B must be in the same country as feature in A in order for it to be joined.  The difficulty, is that the closest feature in B to feature in A may in fact be in a different country, but the next closest feature in B may be in the same country, and this is the feature I'd like to have the spatial join return.  While model builder seems to allow me to iterate through selection of one feature class (the "feature selection" iterator), it does not seem to allow me to iterate a feature selection that draws on attributes for two feature classes, so that I could set the country fields equal before performing the spatial join.  (I'd like to avoid using python since I don't have any experience with it).


texts: Get the distance between 2 WGS84 GEOSGeometries

texts:
I currently want to check if 2 given GEOSGeometry objects of type Point  are within 100m of each other. (length flexible)
Looking at previous asked questions on this topic, most of the answers were totally not specific and even worse sent the author to other frameworks.

My goal is to stay with GEOS / GDAL (geo django) to solve this problem.

The current code I use:

from django.contrib.gis.geos import Point, GEOSGeometry
pnt = GEOSGeometry('SRID=4326;POINT(40.396764 -3.68042)')
pnt2 = GEOSGeometry('SRID=4326;POINT( 48.835797 2.329102  )')
pnt.distance(pnt2) * 100

Found on gis.SO.

It is very difficult for me to confirm that this code works, because I would have to check the distances between to points and check the results against each other.

The main goal is to take 2 points and get the distance in meters based on the WGS84 projection and the GEOS or GDAL framework


texts: Map Package is not showing points on map

texts:
I am fairly new to GIS, but I am working on an interactive legend and I need to just be able to add a map package to the xaml and it place the points on the map. The legend seems to be working. All of the label names are right, but the points do not show up. When I did this before with just each individual layer it worked. Any help would be appreciated. Here is my xaml.

<esri:Map x:Name="MyMap" Background="#FFE3E3E3" WrapAround="True" MouseClick="MyMap_MouseClick" MouseMove="MyMap_MouseMove" MouseRightButtonDown="MyMap_MouseRightButtonDown">
                <esri:ArcGISLocalDynamicMapServiceLayer ID="Calvert_City" Path="C:\Users\jessical\Documents\ArcGIS\CalvertCity_Test.mpk"/>
                <esri:ArcGISTiledMapServiceLayer ID="MyLayer" Url="http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer" />
</esri:Map>

            <Border Background="{StaticResource PanelGradient}" BorderThickness="1" CornerRadius="5"
            HorizontalAlignment="Left"  VerticalAlignment="Top"
            Margin="20" Padding="5" BorderBrush="Black" >
                <esri:Legend Map="{Binding ElementName=MyMap}"

                         LayerItemsMode="Tree"
                         ShowOnlyVisibleLayers="False"
                         Refreshed="Legend_Refreshed">
                    <esri:Legend.MapLayerTemplate>
                        <DataTemplate>
                            <StackPanel Orientation="Horizontal">
                                <CheckBox Content="{Binding Label}"
                                    IsChecked="{Binding IsEnabled, Mode=TwoWay}"
                                    IsEnabled="{Binding IsInScaleRange}" >
                                </CheckBox>
                                <Slider Maximum="1" Value="{Binding Layer.Opacity, Mode=TwoWay}" Width="50" />
                            </StackPanel>
                        </DataTemplate>
                    </esri:Legend.MapLayerTemplate>
                    <esri:Legend.LayerTemplate>
                        <DataTemplate>
                            <CheckBox Content="{Binding Label}"
                            IsChecked="{Binding IsEnabled, Mode=TwoWay}"
                            IsEnabled="{Binding IsInScaleRange}" >
                            </CheckBox>
                        </DataTemplate>
                    </esri:Legend.LayerTemplate>
                </esri:Legend>
            </Border>

If you need anymore information, please let me know. Thanks!


texts: Mapfish 3 variables in scalebar processor

texts:
I would like to be able to change the units of the scalebar processor based on values sent in spec.json - It seems this should be possible from reading the docs, but I can't get it to work.

A snippet of my config.yaml below showing the relevant parts. The config does work until I attempt to use the ${mapUnits} variable instead of a distinct value.

config.yaml

      title: !string
        default: ""
      mapUnits: !string
        default: "m"
      scalebar: !scalebar
        width: 180
        height: 31
        createSubReport: false
        default:
          backgroundColor: "rgba(255, 255, 255, 0)"
          unit: ${mapUnits}
          geodetic: true
          fontSize: 8
          barSize: 8
          padding: 5
          intervals: 2
          align: center
          type: line

I am also declaring 'mapUnits' in the spec:

spec.json snippet

{
    "layout": "A4 Landscape",
    "outputFormat": "pdf",
    "attributes": {
        "mapUnits": "ft",
        "title": "My Map",
        "map": {
            "bbox": [-1910349.8815163244, 6214452.698830085, 1205906.228840141, 8182614.4527394315],
            "dpi": 150,
            "layers": [/*layer config*/]
        }
    }
}

Running this just produces the error Exception in thread "main" java.lang.IllegalArgumentException: invalid unit: ${mapUnits}

Is what I am attempting to do possible?


texts: QGIS combine Vectorlayer and Rasterlayer to one shapefile

texts:
I have to create some Thematic Maps, especially Choropleth Maps on Matlab.

Therefore i've found a shapefile which I can use, but need to draw some lines there.

I've imported my shapefile as a Vector Layer and my .jpeg file as a Raster layer and want to combine both to one shapefile.

I did "Polygonize" my Raster layer and combined both Vector layers to one shapefile, but it doesnt look the way as i want it.
Here are some screenshots:

My Vector Layer:

My Raster layer:

I want this result:

Does anyone know, how I can achieve this result with this transparency?

Best regards.


texts: Getting LiDAR point cloud point spacing using PDAL

texts:
I am trying to use PDAL (in Python) to extract metadata from a LiDAR point cloud (.las).  I am using the "pdal info" command and applying different flags (e.g. -metadata, -stats) to get additional information on the dataset, but I would like to know if it is possible to obtain the point spacing of a LAS point cloud using PDAL?
I initially thought I had the correct information in the scale_x, scale_y, scale_z values returned via "pdal info -metadata", but upon reading more about PDAL my understanding now is that these are scaling factors and not the actual point spacing.
Is point spacing actually retrievable via PDAL?


texts: Prevent loss of timestamp when mapping over image collection

texts:
I am using the below code (taken from this post) to make monthly composites of cloud masked S2 data but the reduce.median won't work because my ImageCollection doesn't have a timestamp. Does anyone know how I can add a timestamp to the masked ImageCollection?

/* Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

var geometry = geometry;

var col = ee.ImageCollection('COPERNICUS/S2')
          .filterBounds(geometry).filterDate('2018-01-01','2019-01-01')
          // Pre-filter to get less cloudy granules.
          .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
          .map(maskS2clouds)

print(col)

//Monthly Images
var months = ee.List.sequence(1,12);
var years = ee.List.sequence(2018,2018);

// monthly composite
var monthlyImages =  ee.ImageCollection.fromImages(
  years.map(function (y) {
  return months.map(function(m){
  var w = col.select('B4', 'B3', 'B2').filter(ee.Filter.calendarRange(y, y, 'year'))
           .filter(ee.Filter.calendarRange(m, m, 'month'))
           //.mean()
           .median()
           //.min()
           .set('system:time_start',ee.Date.fromYMD(y,m,1));
  return w.set('year', y)
           .set('month', m)
           .set('date', ee.Date.fromYMD(y,m,1))
           .set('system:time_start',ee.Date.fromYMD(y,m,1))

});
}).flatten());

print(monthlyImages)

Error: ImageCollection (Error)
reduce.median: Can't apply calendarRange filter to objects without a timestamp.


texts: MapInfo Pro 12.5 Pan Control not going upward when measuring a distance

texts:
MapInfo Pro 12.5  (32 bit)
In one workspace I am measuring a distance - I can pan left right and down, but not upwards.
The pan works in another workspace in all directions.
Why does it not work in the workspace?
Is there a particular setting in the Map Preferences or Map Options that needs to be selected to enable the pan to work in all directions.....?

Below is an image of Map Options


texts: How to compile Mapnik 2.2.0 with mapnik-vector-tile headers

texts:
I want to generate vector tiles with Mapnik 2.2.0 (That's the latest version that can be installed on debian wheezy). According to this answer, I need to compile mapnik with mapnik-vector-tile which provides C++ headers. v0.6.x of mapnik-vector-tile supports Mapnik 2.2.0 so I can use that.

How do I compile mapnik with those headers?


texts: Raster displaying correctly, but incorrect values

texts:
I have a raster dataset created from the pixel classification of a satellite image. The pixels have 4 possible values, representing different land cover classes : 1, 2, 3 and 4. Using ArcMap 10,4, The display works fine, as you can see in the first image.

However, when looking at this raster's properties wanting to change the colors to something "land-cover-ish", the minimal value of the file is not 1, but rather 1,33244 (see image below). I am therefore unable to display the 4 different classes, neither to perform zonal statistics as I intended to do.
Strangely enough, the four classes appear as expected in QGIS.
Has anyone stumbled upon this issue?


texts: “User memory limit exceeded” using Random Forest classification algorithm in Google Earth Engine

texts:
I have see this conversation here
I have exactly the same problem with Random forest classifier and I would like know if it was resolved ? if so, which means did you find to resolve it ?


texts:
Below is a sample one polygon geometry type census shapefile information. Looking to pass a point via (latitude, longitude) so that I can perform a point in polygon and get the details of the polygon using a SQL script.

ID  1
STATEFP 25
COUNTYFP    1
TRACTCE 11200
AFFGEOID    1400000US25001011200
GEOID   25001011200
NAME    112
LSAD    CT
ALAND   13422840
AWATER  681467
geom    0xAD100000010427000000A696ADF545D64440C57


texts: A lot of pointers to load on map-A better solution

texts:
I have to project a lot of points present in my Database say more than 100000. I heard that using geoserver is a better solution as it can render any number of points and from http://swingley.appspot.com/maps/olpts it is almost confirmed that I cannot use the normal code to project many points. I recently heard about cluster strategy that is present in OpenLayers and I think it can help me but I am not sure.

Can anyone help me about the better solution of these two GEOSERVER or CLUSTER strategy and reasons why one is preferred over the other? Remember I have to display many layers on the map.


texts: GeoToolKit conversion from lat/long to UTM

texts: What is the proper way to convert lat/longs to a UTM coordinate system using the Java geotoolkit API? I have not worked with GIS software before, so maybe what I am trying to do is not directly possible. I found a few tutorials, but nothing with the full solution.


texts: Alternatives for the Emerging Hot Spot Analysis/Mann-Kendall statistic (less than 10 time steps)

texts:
I came across ESRI's Emerging Hot Spot Analysis (Space Time Pattern Mining toolbox) and really like what it does. Unfortunately for me, it relies on the Mann-Kendall statistic which apparently needs a minimum of ten time steps to run the analysis. I have less than ten time steps.
I was wondering if anyone knew of an alternative analysis that can be used with fewer than ten time steps?
I can execute the analysis in Arc or in R.
Many thanks in advance.


texts: How does ArcGIS Desktop match WMS EPSG codes with ESRI CRS names and parameters and how to add a missing one?

texts:
In ArcGIS Desktop 9.3.1, when connecting to a WMS layer, the server sends to ArcGIS the supported CRSs (coordinate reference systems) as a list of EPSG codes.

I wonder how ArcGIS (ArcMap and ArcCatalog) recognize the coordinate reference system name and related parameters for each EPSG code in the list sent by the server, because I need to add one missing EPSG code -> coordinate reference system correlation in ArcGIS in order to workaround a bug that occurs in the following scenario:

when the server returns an EPSG codes list in which the first EPSG code is not recognized by (does not match any coordinate system name known by) ArcGIS

[EPSG list returned by the server in the GetCapabilities command replay]

<CRS>EPSG:6706</CRS>
<CRS>EPSG:4258</CRS>
<CRS>EPSG:25832</CRS>
<CRS>EPSG:25833</CRS>
<CRS>EPSG:25834</CRS>

[Layer Properties / Source / Data Source excerpt from ArcMap]

Supported Coordinate System(s):
 EPSG:6706   -
 EPSG:4258   - GCS_ETRS_1989
 EPSG:25832  - ETRS_1989_UTM_Zone_32N
 EPSG:25833  - ETRS_1989_UTM_Zone_33N
 EPSG:25834  - ETRS_1989_UTM_Zone_34N

then ArcMap fails to draw the layer, even changing the data frame coordinate system in one of the others (both from the Data Frame Properties or from the Change Coordinate System dialog).


texts:
I am exploring ArcGIS Desktop 10.6 as I am mostly used to QGIS. I am making use of the world light grey canvas as this is the most suitable for me.

However, I would prefer if the oceans appear blue while every other feature is grey. I can't seem to find where I can make this kind of change. The topographic one has too much information especially when you zoom into it.

How can I make these edits?


texts: python queries distance the nearest object

texts:
I have 2 point layers, I want to calculate the shortest distance between 2 point layers by python code and enclose the conditions of satisfaction.
For example, a point a when there is the shortest distance with point b and this a point must have  a field that satisfies the given condition.

Is there any python to do it?


texts: undefined reference to `TIFFDefaultStripSize@LIBTIFF_4.0'

texts:
I am trying for the first time to use gdal in Ubuntu 15.10 with gcc. I call gdal.h in my program.
At first I got the error "undefined reference to `GDALOpen'", so I check online and I found that I had to link gdal library.
So, now the command line is

gcc prova.c -lm -I/usr/include/gdal/ -L/usr/lib -lgdal

And the error I get this time is

/usr/lib/x86_64-linux-gnu/libpoppler.so.52: undefined reference to TIFFDefaultStripSize@LIBTIFF_4.0
/usr/lib/x86_64-linux-gnu/libpoppler.so.52: undefined reference to TIFFFdOpen@LIBTIFF_4.0
/usr/lib/x86_64-linux-gnu/libpoppler.so.52: undefined reference to TIFFWriteScanline@LIBTIFF_4.0
/usr/lib/x86_64-linux-gnu/libpoppler.so.52: undefined reference to TIFFClose@LIBTIFF_4.0
/usr/lib/x86_64-linux-gnu/libpoppler.so.52: undefined reference to TIFFSetField@LIBTIFF_4.0

I am really a newbie in C, I am still learning.
I was not able to find any solution online, could anyone help me?


texts: Cloning a transactional version in ArcGIS

texts:
In ArcGIS documentation, I cannot find any references to cloning a transactional version. Is there a built in method for it?

I think I know how to do it manually in Arcpy, basically creating a new state branching from a parent and copying all the deltas with the given state (clone source) to the new state (clone) and of course adding the new state to the version tree under the parent, but I am wondering if there is a more neat way of doing it.

Update: Lets put it in more abstract way: Assume version hierarchy like A->B, A->C and C->D. i.e. B and C are decedents of A and D is decedent of C. The objective is to get an independent copy of C and call it E, while keeping the proper hierarchy. That is, after cloning, we will have A->E and E will contain exactly all the changes in C.
Note that, creating a version off of C will not cut it b/c if C changes, those changes will be reflected in E, which we do not want. Also, taking a branch off of A will not work since it does not have what is in C. Of course, we do not want to reconcile C into A and then get a version off of A - that is against the whole goal of having all these versions around.


texts: My QGIS project won't open - unable to unzip .qgz file

texts: I created a QGIS project yesterday, and saved it in one of my user's folder. The file extension is qgz. Today I can't open it anymore. In the project home folder in QGIS it says Unable to unzip file and then file path where the file is.


texts:
I'm completely new to GDAL and I'm going to use it to work with satellite images, DEMs, etc.

I have some experience working with Visual C++ and the OpenCV library. And in fact my goal is to combine OpenCV, GDAL and OpenGL for developing software for remote-sensing purposes. The programming language will be Visual C++.

My problem is that in searching the net, I couldn't find a good book for learning GDAL and in fact most of the open source codes found in the Net are written in the Python language but I'm not familiar with this language and I don't want to get familiar.  Surely I know that the best language for programming with GDAL is Python but the rules of the company where I work force me to use Visual C++.

So with the conditions mentioned above, what do you suggest as the best way to study GDAL?

The first piece of program that I want to write would be:

Read a satellite image file with the help of the GDAL library.
Extract image's header information.
Write a program to  apply a Gaussian filter on the image that was read in step 1.
Convert the image read in step 1 to an image of the OpenCV family to do some other processing.


texts: Accessing PostGIS layers as non-data owner via QGIS?

texts:
I'm trying to create a database user (postgres) to access PostGIS data via QGIS.

I've applied the following commands:

GRANT ALL ON SCHEMA PUBLIC TO ene;
GRANT SELECT ON ALL TABLES IN SCHEMA public to ene;
GRANT ALL PRIVILEGES ON DATABASE "raiateaBD" TO ene;

The access tab for this database looks like this:

However, when creating a connection for this user in QGIS, the tables are not available.

Logging in as postgres user, the tables are visible and available to be loaded:

Furthermore, a project previously created with data loaded by postgres user (data owner) provokes the handle bad layers window when logged in as ene user.

The QGIS log has this error for each of the layers that are in the database:

I added the following, and the user 'ene' still does not have access:

GRANT ALL ON TABLE public.spatial_ref_sys TO ene;
GRANT ALL ON TABLE public.geography_columns TO ene;
GRANT ALL ON TABLE public.geometry_columns TO ene;

BREAKTHROUGH

Layer is accessible if I call:

GRANT SELECT ON TABLE public."table_name" TO ene;

However, I'd like to not calling it one by one, but the following command is not doing it:

GRANT SELECT ON ALL TABLES IN SCHEMA public TO ene;

What am I missing?


texts: How to run a tool from Whitebox GAT in ArcGIS?

texts:
I have a series of tools that I need to run in ArcGIS to hydrologically condition my DEM, but the first step in this process is running the Breach Depressions tool from Whitebox GAT (I have found that this tool works really well for breaching streams across bridges and dams).

Is it possible to work entirely in ArcMap and call the breach depression tool from whitebox? I tried executing the python code in ArcGIS but it doesn't seem to work.
Here is what I tried:

wd = "D:\LIDAR\02060006_Patuxent\01_HeadwatersPatuxent\0104_HUC12"
pluginHost.setWorkingDirectory(wd)
inputFile = wd + "DEM_clip.dep"
outputFile = wd + "DEMBreach.dep"
maxBreachLength = "50"
outputPointer = "true"
args = [inputFile, outputFile]
pluginHost.runPlugin("BreachDepressions", args, False)

I'm pretty new to python and coding in general so the answer may be more obvious than I think. I should add that right now I'm running my tools in a geoprocessing model so ideally I would like to somehow add the whitebox breach depressions tool into my geoprocessing model, but would be open to trying it out in a python script.


texts:
I use Leaflet to show a WMS Layer via a proxy.php (for authentication). That is working on most WMS Layers. But I have some WMS Sources, which are having a small shift in the tiles. See Overlay Map in the example.

<!DOCTYPE html>

    TEST MAP
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0">

<link rel="stylesheet" href="https://unpkg.com/leaflet@1.6.0/dist/leaflet.css" />
<script src="https://unpkg.com/leaflet@1.6.0/dist/leaflet.js"></script>

var map1  = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 30,
    maxNativeZoom: 18
});

var map4 = new L.tileLayer.wms("proxy.php?url=http://geoshop.kommunalportal.ch/SWMS?", {
    layers: 'lk_wms_wasser',
    format: 'image/png',
    transparent: true,
    maxZoom:30
});

var map = L.map('mapid').setView([47.549904, 9.152927], 18);

var baseMaps = {
    "map1": map1

};

var overlayMaps = {
    "map4": map4
};

L.control.layers(baseMaps, overlayMaps).addTo(map);

map1.addTo(map);
map4.addTo(map);


texts: Creating Variable Distance Buffers using QGIS

texts:
I would like to create variable distance buffers for some points in QGIS 3.4.

I have seen in previous posts that there used to be a tool for this in older versions of QGIS, but that this tool wasn't continued in the following update.

Is there a way to create buffers based on information from a column of the attribute table?


texts: Optimize sampling on band 1 of a big dataset

texts:
We do MANY samplings on band 1 of a big dataset, using rasterio dataset.sample function. (https://github.com/mapbox/rasterio/blob/master/rasterio/sample.py)
Our python code looks something like this:
xy_list = [(448758.24921202223, 6239730.062386721)] # normally many more
with rasterio.open(tif_file_path) as dataset:
  sample_result_tuple = list(dataset.sample(xy = xy_list , indexes = 1))
  ...

Is there any way to speed this up, maybe by using some sort of fancy indexing or preprocessing? Any suggestion at all is welcome, even if it is not in python.


texts:
I work with raster layers (tif). I have to seperate buildings and vegetation (trees) from a DSM.

DSM - DTM with the rastercalculater gives me the buildings and the trees. (DSM and DOM combined)

But then I get where no trees or buildings are numbers like -0.00332... how can I "delete" those values from my new raster that all values under 2 are = 0?


texts: Why does using next() with an update cursor nested inside a Search cursor yield such weird results

texts:
I did this:

with arcpy.da.SearchCursor(
        in_table = in_table, field_names = '*', where_clause = "",
        sql_clause = (None, "ORDER BY date"))\
        as ordered_cursor:

        with arcpy.da.UpdateCursor(
            in_table = in_table, field_names = '*')\
            as cursor:
            for row in cursor:
                row = next(ordered_cursor)
                cursor.updateRow(row)

I expected this to sort my table by date. Something like this,

1 3 5 1/02/1300
6 9 7 1/03/1300
9 9 8 6/09/1342
7 7 7 7/07/1777

What I got was something along the lines of,

1 3 5 1/02/1300
1 3 5 1/02/1300
6 9 7 1/03/1300
1 3 5 1/02/1300

Why did this happen? Is ordered_cursor still iterating over the values of in_table? I noticed that I was unable to delete the table after I did this while I had the python window that the code executed in open. This was the result of a schema lock. I just want to clarify that this code executes without raising an exception.

I tried using sort management and it doesn't work where I'm using Enterprise connections. I'm not looking for a way to make sort management work. I'm also not looking for a way to sort things. I'm just looking for an explanation of this behavior.


texts: Draw shape in map and save to server ArcGIS JavaScript version 4+

texts: is it possible to use sketchViewModel to draw polygons and using applyEdits() to save as feature layer in server?


texts: How to select only the lines that has points at at least one end?

texts: I have 2 separate feature classes, line and point. I want to select line feature class based on the location of points on a point feature class. If the points are on the line itself then I don't want to select the line. But, if they are at its end/start, the line should be selected.


texts: Openlayers Transform not transforming a point from EPSG 27700 to EPSG 4326

texts:
I am trying to add Google street view functionality to a Openlayers Map. I have started by using 3 divs to hold the Openlayers Map (which uses a local tilecached base map), a google api map and another for street view.

If I try to transform the lonlat from EPSG:27700 to EPSG:4326 so that they can be passed to the google map object as the center location, the coordinates are not transformed.

I have the following scripts included in the html:

    <script src="script/openlayers/proj4js-compressed.js"></script>
    <script src="script/openlayers/OpenLayers.js"></script>
    <script src="script/openlayers/defs/EPSG27700.js"></script>
    <script src="http://ajax.aspnetcdn.com/ajax/jQuery/jquery-1.8.1.min.js"></script>
    <script src="https://maps.googleapis.com/maps/api/js?v=3.exp&sensor=false"></script>
   <script src="script/map.js"></script>

If I alert the lonlat variable after transforming it I get the original EPSG 27700 value of lon = 311711, lat = 530372, rather than the converted EPSG 4326 value.

The Openlayers map is set up here:

    //private variables
    var map = null;
    var streetview = null;
    var osgrid = new OpenLayers.Projection("EPSG:27700");
    var googlegrid = new OpenLayers.Projection("GOOGLE");
    var mapextent = new OpenLayers.Bounds(295000.0,505000.0,340000,565000.0);

    var options = {
         projection: osgrid,
    displayProjection: osgrid,
    maxExtent: mapextent,
    restrictedExtent: mapextent,
    units: 'm'
    }
    var osoptions = {
    scales: [
        250000,
        107716.5936,
        70866.18,
        47243.9527,
        22677.1776,
        15000.236,
        10000.9708,
        6000.309,
        3000.91484,
        2000,
        500
    ]
};

    //allow testing of specific renderers via "?renderer=Canvas", etc
    var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
    renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;

    //Set mapping this object for use in internal functions
    var that = this;

    var addControls = function() {
           var nav = new OpenLayers.Control.Navigation();
           var zoomBox = new OpenLayers.Control.ZoomBox();
           var layerSwitch = new OpenLayers.Control.LayerSwitcher();
           map.addControls([nav, zoomBox, layerSwitch]);
    };

    //Public Objects
    this.tmsurl = ["http://servername1/tilecache/",
                    "http://servername2/tilecache/",
                    "http://servername3/tilecache/",
                    "http://servername4/tilecache/"
                ];

    this.popup = null;

    this.init = function() {

        map = new OpenLayers.Map("map_canvas", options);

        var osmap = new OpenLayers.Layer.TileCache(
            "OS Map",
             that.tmsurl,
            'base_StandardMap',
            {
                 options: osoptions,
                format: 'image/jpeg',
                isBaseLayer: true
            }
        );

        addControls();

         map.addLayers([aerialmap, osmap]);

        var lon = 311711;
        var lat = 530372;
        var zoom = 0;

        map.zoomToExtent(mapextent);

        that.centreMap(lon, lat, zoom);

The Google Map and Street view is set up here

var lonlat = new OpenLayers.LonLat(lon, lat);
lonlat.transform(osgrid,new OpenLayers.Projection("EPSG:4326"));
        var location = new google.maps.LatLng(lonlat.lat, lonlat.lon);
        var mapOptions = {
          center: location,
          zoom: 14,
          mapTypeId: google.maps.MapTypeId.ROADMAP
        };
        var gmap = new google.maps.Map(
            document.getElementById('gmap'), mapOptions);
        var panoramaOptions = {
          position: location,
          pov: {
            heading: 34,
            pitch: 10,
            zoom: 1
          }
        };
        var panorama = new  google.maps.StreetViewPanorama(document.getElementById('streetview'),panoramaOptions);
        gmap.setStreetView(panorama);

    },


texts: Using ST_ClusterDBSCAN on results of ST_ClusterKMeans for nested clustering with PostGIS?

texts:
I would like to pipe the results of ST_ClusterKMeans() into a ST_ClusterDBSCAN() query.

The ST_ClusterKMeans gives me a fairly good results as shown in this screenshot. I classified the union of the point geometries returned from the grouped cluster and I overlaid the ST_PointOnSurface() geometries from the multipoints as black circle (size is dependent on the number of clustered points).

This is the query which returns the multipoint geometries:

SELECT
      ST_Union(cluster.geomcntr) AS geom,
      count(cluster.geomcntr) AS c
    FROM (
      SELECT
        qid AS id,
        brand AS brand,
        store AS label,
        --ST_ClusterKMeans(geomcntr, 200) OVER () AS cid,
        ST_ClusterDBSCAN(geomcntr, 0.1, 1) OVER () AS cid,
        geomcntr
      FROM retailpoints
      WHERE
        ST_DWithin(
          ST_MakeEnvelope(
            -17.951660156250004,
            59.512029386502704,
            9.953613281250002,
            49.439556958940855,
            4326),
          geomcntr,
          0.00001
        )
    ) cluster
    GROUP BY cluster.cid;

I would like to nest this query inside a query that calculates the ST_ClusterDBSCAN cluster inside each of the multigeometries, effectively splitting KMEANS cluster which are very sparse into separate village cluster.

Edit: Dan Baston's suggestion works a charm. Does exactly what I was trying to do. Here is a screenshot with the PointOnSurface centroids from the DBSCAN cluster which use the KMEANS cluster as input.


texts: Finding a more efficient way to execute a long-run python script

texts:
I'm executing a python script which does these steps:

iterates through 19 features (makes pairs "each one with each other", order doesn't matter: A-B = B-A)
Creates Closest Facility Layer
Adds location points
Solves Closest Facility route
Exports into single .shp for each pair

I need to get a polyline for each of the pairs of 19 combinations. This will result (if I count correctly) in 19*18/2 = 174 shapefiles.
When I run the script, it works fine until about 40th pair. Then it gets slower and slower. I let it run for about 3 hours and I got 142th pair...then ArcMap crashed with an error screen saying all my work is lost. If it ended alright I would then merge all the shapefiles into a single one with 174 features in it.

Is there any way to make the script more efficient or a way how to solve this problem?

This is the script:

 from collections import defaultdict
 import itertools
 import pprint

 import arcpy
 from arcpy import env
 import os

 ddict = defaultdict(list)
 for start,end in itertools.combinations(range(19), 2):
ddict[start].append(end)

try:
arcpy.CheckOutExtension("Network")
env.workspace = r"D:/Dokumenty/GISdata/PID"
env.overwriteOutput = True

inNetworkDataset = r"D:/Dokumenty/GISdata/PID/New/silnice_useky_krnap_ND.nd"
impedanceAttribute = "Length"
accumulateAttributeName = ["Length"]

for key,value in ddict.items():
    for v in value:
        inFacilities = r"PID/New/Obce_body_luzka" # melo by byt key
        inIncidents = r"PID/New/Obce_body_luzka" # melo by byt v
        outNALayerName = "Dvojice_bodu_{}_{}".format(key, v)
        outLayerFile = os.path.join(r"D:/Dokumenty/GISdata/output2", outNALayerName + ".lyr")
        NAResultObject = arcpy.na.MakeClosestFacilityLayer(inNetworkDataset,outNALayerName,impedanceAttribute,"TRAVEL_TO","",1, accumulateAttributeName,"NO_UTURNS")

        outNALayer = NAResultObject.getOutput(0)

        subLayerNames = arcpy.na.GetNAClassNames(outNALayer)
        #Stores the layer names that we will use later
        facilitiesLayerName = subLayerNames["Facilities"]
        incidentsLayerName = subLayerNames["Incidents"]

        import arcpy

        intable = "D:\Dokumenty\GISdata\PID\New\Obce_body_luzka.dbf"

        fields = arcpy.ListFields(intable)

        fieldinfo = arcpy.FieldInfo()

        for field in fields:
            fieldinfo.addField(field.name, field.name, "VISIBLE", "")
        arcpy.MakeTableView_management(intable, "temp_table", "FID = {}".format (key), "", fieldinfo)
        arcpy.MakeTableView_management(intable, "temp_table2", "FID = {}".format (v), "", fieldinfo)
        fieldMappings = arcpy.na.NAClassFieldMappings(outNALayer, facilitiesLayerName, False, fields)
        arcpy.na.AddLocations(outNALayer, facilitiesLayerName, "temp_table", fieldMappings, "")
        arcpy.na.AddLocations(outNALayer, incidentsLayerName, "temp_table2", "", "")

        arcpy.na.Solve(outNALayer)

                arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"RELATIVE")
        routesLayer = arcpy.mapping.ListLayers(outNALayer, "routes")[0]
        arcpy.FeatureClassToShapefile_conversion(routesLayer, r"D:/dokumenty/gisdata/output4")

print "Script completed successfully"

except Exception as e:

import traceback, sys
tb = sys.exc_info()[2]
print "An error occured on line %i" % tb.tb_lineno
print str(e)


texts: ArcGIS/Server can not display the digital places as ArcGIS/desktop

texts:
The number of decimal places can be set in ArcGIS desktop as shown in the screenshot.

However, ArcGIS server fails to do this.

What could be the issue here? How can I reflect the displayed digital places from the desktop to the server?


texts:
I'm having difficulty troubleshooting why a number of my mapfile layers aren't rendering. I'm a GIS noob, so I'm probably missing something obvious and apologize in advance; I just don't know where else to turn to.

I built a map project in QGIS and arranged the layers how I'd like; everything displays fine there. But when I attempt to query MapServer for those same layers, only some of them are rendering.

To illustrate, here's an example of how I'm using the mapfile:

MAP
NAME "WASHINGTON_COUNTY"

SIZE 800 600
EXTENT 721132.369589 1497660.070845 825840.827167 1814625.359158
UNITS METERS
SHAPEPATH "geo"
TEMPLATEPATTERN "washington"
IMAGECOLOR 255 255 255
IMAGETYPE PNG

WEB
    HEADER templates/header.html
    TEMPLATE "set in index.html"
    FOOTER templates/footer.html
    MINSCALE 1000
    MAXSCALE 2550000

    IMAGEPATH "C:/ms4w/tmp/ms_tmp/"
    IMAGEURL "/ms_tmp/"
END

LAYER
    NAME tiger_counties
    DATA "tiger_counties"
    STATUS OFF
    TYPE POLYGON

    CLASS
        NAME "county"

        STYLE
            COLOR 200 200 200
            OUTLINECOLOR 222 222 222
        END
    END
END

LAYER
    NAME ny_washington_centerlines
    DATA "ny/washington/centerlines"
    STATUS OFF
    TYPE LINE

    CLASS
        NAME "road"

        STYLE
            COLOR 255 255 0
            OUTLINEWIDTH 1
            OUTLINECOLOR 0 0 0
        END
    END
END
END

And here's the image rendered from the query string:

"map_size=592+886&img.x=441&img.y=257&mode=browse&layer=ny_washington_centerlines&layer=parcels&zoomdir=0&zoomsize=2&imgxy=400.0+300.0&imgext=718662.996883+1676819.127354+824362.189500+1756060.449478&map=C%3A%2Fms4w%2Fapps%2Fwashington_county_mapserver%2Fwashington_county.map&root=%2Fwashington_county_mapserver%2F&map_web=+TEMPLATE+washington_main.html"

:

Now here's a map of the same region, as rendered in QGIS:

What is QGIS doing by default that MapServer doesn't?


texts:
I want to create oval buffers 450 by 350 m tilted by 45° using the v.buffer-tool from the GUI.
When I enter everything necessary and click run, in return I get the following error message (same pops up for other Grass-algorithms):

2019-01-29T11:56:21     CRITICAL    Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 415, in processAlgorithm Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 372, in executeGrass for line in iter(proc.stdout.readline, ''):
File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\codecs.py", line 322, in decode (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xf6 in position 81: invalid start byte

A possible reason could be the following warning I get upon startup of QGis (it doesn't pop up, but it is written in the protocol):

2019-01-29T11:43:18     WARNING    Konnte
C:/PROGRA~1/QGIS3~1.4/apps/qgis/plugins/grassplugin7.dll nicht laden
(Grund: Die Bibliothek
C:\PROGRA~1\QGIS3~1.4\apps\qgis\plugins\grassplugin7.dll kann nicht
geladen werden: Das angegebene Modul wurde nicht gefunden.)
2019-01-29T11:43:18     WARNING    Konnte
C:/PROGRA~1/QGIS3~1.4/apps/qgis/plugins/grassprovider7.dll nicht laden
(Grund: Die Bibliothek
C:\PROGRA~1\QGIS3~1.4\apps\qgis\plugins\grassprovider7.dll kann nicht
geladen werden: Das angegebene Modul wurde nicht gefunden.)
2019-01-29T11:43:18     WARNING    Konnte
C:/PROGRA~1/QGIS3~1.4/apps/qgis/plugins/grassrasterprovider7.dll nicht
laden (Grund: Die Bibliothek
C:\PROGRA~1\QGIS3~1.4\apps\qgis\plugins\grassrasterprovider7.dll kann
nicht geladen werden: Das angegebene Modul wurde nicht gefunden.)

Basically it says "GRASS wasn't installed". So, how do I solve this? Just reinstall QGis?


texts: Which geographic coordinate system should I use to work in meters in Angola?

texts:
I am using WGS84 in Qgis 2.18, but when I try to use some Qgis Plugins (like Profile Tool and Q-Ras) they don´t work. Is it a coordinate problem? Or something else?

Profile Tool:

It makes the profile but the horizontal scale instead of being in meters is something strange (0,002...0,004...etc). The vertical scale is ok.

Q-Raz:

Error message states:

  2018-01-08T18:04:25   1   Traceback (most recent call last): File
  "C:/Users/Jorge Matos/.qgis2/python/plugins\Qgis2Ras\qras_dialog.py",
  line 113, in runProfile
  main(self.river,self.XSection,self.textfile,self.dem,resolution)
  File "C:/Users/Jorge Matos/.qgis2/python/plugins\Qgis2Ras\RASout.py",
  line 284, in main outfile = open(file_strin_path, 'w') IOError: [Errno
  22] invalid mode ('w') or filename:


texts: Quering a ArcSDE featureclass attribute table after joining

texts:
I am trying to query a ArcSDE Feature Class after joining using Iqueryfilter

I have two SDE feature classes.
1)  Subdivisions :      which  has details of subdivision and districts (parent)

2)  Districts:   Which has details of districts and States (parent)

I want to find out what are the subdivisions in a state. For this I joined Districts table to Subdivision table and  now I can do attribute query in ArcMap and get the divisions in a state using select by attribute tool (example :  Districts.statename =  'A'  selects all the divisions in the state ‘A’)

Now I want to do the same using VBA -   I am using the  following method but it gives me  ORA-00904 error :  “Districts”.”statename” : invalid identifier.

Anybody have faced same type of issue erliar?  What is best way to get this done.
Thank you in advance.

Dim pQFilter As IQueryFilter

Set pQFilter = New QueryFilter

Dim pFeature As IFeature

pQFilter.WhereClause =  " Districts.statename = 'A'"

Dim myFeatureLayer As IFeatureLayer

Set myFeatureLayer =   “Districts”

Dim pFCursor As IFeatureCursor

Set pFCursor = myFeatureLayer.Search(pQFilter, False)

Set pFeature = pFCursor.NextFeature()


texts: How to do diagram block with open source software?

texts: I would like to do a block diagram (topographic and geological representation). Is it possible to do this with QGIS (or other open source software), in Windows 7?


texts: Copying vertex coordinates from one line into another with same geometry type

texts:
I have the following QGIS version.

How is possible to copy the vertex data information from a line into another line? Maybe a certain plugin shall be used for this purpose.

The data are 407 points with (X,Y) coordinates tuples. I need to copy this information to another line with the same geometry type.


texts: Building Python script that exports to several output files using ArcGIS Desktop?

texts:
I'm trying to build a Python Script that will receive one feature layer that will be clipped by several adjacent grid squares (each in a different shapefile) and output several shapefiles.

Now, I've built in the past a loop script that receives a list of files and does a single action on each file, like this:

for file in os.listdir(ws):
    if file.endswith(".shp") or file.endswith(".dxf") or file.startswith("str"):
    print file

But I'm not succeeding in making a loop that will receive the input shapefile and clip features (shapefiles) and EXPORT to different shapefiles (I mean how do I change the output shapefile's name for every grid square ?). Can someone show how it is done?

I'm using python 2.6 and Arcgis 9.3.1


texts: Adding Virgin Islands and Puerto Rico to CartoDB?

texts:
I am new to Cartodb.

I have started with the ne_50m_admin_1_states data.

How can I add Virgin Islands and Puerto Rico?

Is there geometry I can get?

Is there an easy way to add them?


texts: Get number of rules in Geoserver style using webservice

texts: How can I know the default style for a layer in Geoserver, and the number of rules it contains, using a REST call or service?


texts:
How do you properly use ST_Retile https://postgis.net/docs/RT_Retile.html

there are no directions for it

tried:

update slopes.dem_slope set rast=st_retile(rast,100,100)

and

select st_retile('slopes.dem_slope','rast',100,100)

and

select st_retile('slopes.dem_slope'::regclass,'rast',100,100)

all throwing errors


texts: st_difference creates two linestrings, how can I keep only one?

texts:
I am using R and the sf package to manipulate shapefiles.
I have transect lines that intersect an island going a few hundred meters inland (see image below). The objective of my project is to trim each transect such that the inland portion is removed.  I used sf::st_difference to clip these transect lines with the island shapefile: Trimmed_Transects<-st_difference(st_geometry(Transects), st_union(st_geometry(Coasts)))
This worked great (see image below)!
However, for rivers, lagoons, and other locations the lines are cut multiple times. For my project, I only want to retain the outermost lines that intersect the coastline (below in red ) and want to remove the inland lines (below in blue). Could anyone recommend a good method to select only the outermost lines?