texts: Importance of WGS84 realization for centimeter-level accurate UTM transformation

texts:
I am working with a high-end GNSS receiver connected to an RTK correction data service. I need to measure the position of a ground point with centimeter-level (2-5cm) accuracy. The GNSS receiver provides the data with respect to the WGS84 datum (as longitude, latitude), which I want to convert to UTM coordinates in the 30N zone (as northing, easting).
I read here https://confluence.qps.nl/qinsy/latest/en/world-geodetic-system-1984-wgs84-182618391.html
that WGS84 has multiple realizations, and I wonder if this has to be taken into account to obtain a centimeter-level accurate position of the point in UTM.
More precisely, my questions are:

Are the outputs of my GNSS receiver automatically referring to the latest WGS84 realization (G1762), or do I have to ask the sensor manufacturer which realization is being used?
Does the transformation from WGS84 to UTM depend on the WGS84 realization? Or is the transformation the same for all realizations and I need to specify the realization used when reporting the ground point's UTM position? If the transformation depends on the realization, how can I specify it e.g. in pyproj or what is an appropriate program to use?


texts: Open Jupyter Notebook with Open Data Cube in docker

texts: I run open data cube (ODC) docker (opendatacube-core and jupyterhub) and I expected that I can run jupyter notebook into it with ODC. However, I didn't. How can I run container (docker) of ODC and I can use Jupyter notebook with ODC?


texts: Processing SQLite files using ArcGIS Desktop?

texts:
I'm using qgis and grass for my work. I've produced a shapefile which is over 100mbs (s a result of some statistical analysis).

If i convert it to sqlite format (spatially enabled) the data shrinks down to 15 mbs. Which is acceptable for sharing it with email or dropbox.

BUT my colleague is using arcgis 9.3 - will he be able to read the data without any problems?


texts: Calling QgsProjectionSelector from QGIS Python plugin?

texts:
I am trying to create a dialog that allows the user to select a coordinate system from either the recently used projections or a general list of all available projections. This would basically offer the same functionality as the standard dialog that pops up when you load in a shapefile with no known CRS.

Is there a way to call the "Coordinate Reference System Selector"-dialog from within a plugin? I know it is possible to query the recently used projections and all available projections, but it seems like a lot of overhead to recreate something that is already there in QGIS.


texts:
I need to extract data from the polygon that is located in the coordinates (lat-long).
I am using QGIS 3.16 for creating projects and QGIS Server for requesting info. I need to extract data included in the layers that are located in a specific coordinates.
For example:

The pink layer is the shapefile and the layer with numbers is WMS.
And I want to know the info for these layers (WMS and shapefile) in the coordinates of the yellow point.
Is there a solution?
UPDATE #1:
The idea is getting the layers' info from an API with an HTTP Request.

With QGIS requests like ...?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&... it can be possible to obtain this information?


texts: Does zonal statistics do weighted average?

texts:
When I compute mean using zonal average in ArcGIS, is it taking into account the 'COUNT' in the value attribute table? e.g. if I have a raster with

VALUE     COUNT
1           10
2           20
3           30

If I run zonal statistics on above raster, will it ignore the COUNT column? Any references in your answer will be helpful


texts: Can't join multiple DEMs using QGIS 'Build Virtual Raster' if they span 2 CRS's

texts:
I have 20 odd DEMs downloaded from https://search.asf.alaska.edu/ (ALOS PALSAR)

Location is the Falkland Islands.

Here is an example just adding 2 of the tiles into QGIS which differ in terms of their CRS

EPSG:32720 - WGS 84 / UTM zone 20S

EPSG:32721 - WGS 84 / UTM zone 21S

As soon as I try putting them into a Virtual Raster using following settings:

I get only a single DEM (20S CRS) within the Virtual Raster which would probably make sense:

Switching Build Virtual raster settings to 21S for 'override projection'

And the same 20S raster just shifts to the right, no sign of 21S tile:

If I turn on 'allow projection difference:

I get the 21S tile flipping the other side of the 20S tile

Is it just not possible to have DEM tiles that cross 2 CRS's within a Virtual Raster, and if not is there another way?


texts: How to edit an existing layer using leaflet?

texts:
I am used to working

with OL3 and I wanted to try using leaflet and see if it is as powerful as OL. I loaded a map with some data from a GeoServer WFS layer using this tutorial from Georepublic.

var group = new L.featureGroup().addTo(map);
var geojsonlayer;
function handleJson(data) {
    //      console.log(data);
    geojsonlayer=L.geoJson(data, {
        style: function (feature) {
            return {
                "weight": 2,
                "opacity": 0.65
            };
        },
        onEachFeature: function (feature, my_Layer) {
            my_Layer.bindPopup("ID : "+feature.properties.id+"<br />Name : "+feature.properties.name);
        }
    }).addTo(group);
    map.fitBounds(group.getBounds());
}

I wanted to add some editing and found the Leaflet.draw plugin I loaded it, and so far I can add a feature to the database ( knowing how painful it will be with WFS-T I used an AJAX call to a php script instead ). Here is how I call the draw control :

var drawnItems = new L.FeatureGroup();
map.addLayer(drawnItems);
var drawControl = new L.Control.Draw({
    edit: {
            featureGroup: drawnItems
    }
});
map.on('draw:created', function (e) {
    //Here I call the php script that stores the created feature in the DB
});

The problem is when I want to edit the existing features using the Leaflet.draw plugin, so I thought of changing the featureGroup in edit to group like this :

var drawControl = new L.Control.Draw({
    edit: {
            featureGroup: group
    }
});

but now I get this error:

Uncaught TypeError: Cannot read property 'enable' of undefined

I have done some research and found that I need to enable editing of the layer with this : layer.editing.enable(); the problem is where to add it so that when I click on the edit icon in the Leaflet.draw control I can edit my layer.

I kind of liked the Leaflet.draw plugin menu on the side, is it possible to solve my problem using this plugin or should I make my own editing menu?


texts: Creating Polygons of Equal Density from Points in ArcGIS Desktop?

texts:
I have a table with latitude and longitude for several thousand points and need to create polygons of equal density around them. Thiessen Polygons seemed like a good route to go, but they only surround a single point.

Another method I found was to create a Near table as shown here, but there are potentially 50,000 or more points to be analyzed, and my computer simply can't handle that.


texts: Statistical comparison between different rasters using R?

texts:
I need to compare, statistically, different rasters with a pattern raster (in this case a Revised Universal Soil Loss Equation (RUSLE) raster). I used different methods to calculate the environment vulnerability, and now, I need to calculate which method is closer to the RUSLE raster. I thought about subtracting each method raster to the RUSLE raster and then compare the means of the resulting rasters.

Anyone have any idea using R?


texts: Google Earth Engine : define polygon Labels for assets

texts:
I'm using the importation assets option to import my shapefile into GEE.
Everything is fine, but I would like to modify the labels of my regions.
My graphs don't have the proper label, and for a table exported with more than 50 regions I don't to do it manually.
For a ee.Geometry.Rectangle(-122.45, 37.74, -122.4, 37.8), {label: 'City'})it is straightforward. But I don't know how to do it for a shapefile.

Here is a description of my shapefile once imported:

type:
FeatureCollection
id:
************/*****************/*************
version:
1518788242243348
columns:
Object (4 properties)
NAME_0:
String
NAME_1:
String
NAME_2:
String
system:index:
String
properties:
Object (1 property)
system:asset_size:
301487


texts: Old project on new PC: bad allocation when opening

texts:
Just got a new laptop (windows7). When opening a QGIS project (shape layers on server), I get an error window "bad allocation.
Project was made in Wroclaw, new PC has Lisboa.

Any suggestions?


texts: Calculating Field in ArcGIS Pro ModelBuilder

texts:
Trying to setup my first automated mapping project and I need to add a field and then set it to 0. I've tried to use the Calculate Value tool, but I don't have any idea how to actually use this as it is completely different than the field calculator. (I will note I've searched extensively for it, but any example I find are more complex than I need and nothing seems to work).

What does the expression and code block need to be to set the entire Column to 0?


texts: Discrepancy between raster calculator and geomorphic change detection elevation change results

texts:
I am using two rasters to determine the elevation change in a stream channel.

When I used the Geomorphic Change Detection Software add-in (version 7.4.1), I got a high elevation change of 8.8 ft and a low of -5.5 ft. I  used the raster calculator to determine the difference between the same two raster layers, I got a high elevation change of 4 ft and a low of -4 ft.

What can be causing such a large difference in these results? Is one more trustworthy than the other?


texts: Remote sensing: How to start with interferometry?

texts:
I need to start learning about the interferometry techniques available/used for remote sensing data (radar or optical).

Which are the main applications?
Which techniques are applied?
Which Earth Observation data is more adequate?
Which software (free or comercial) is better?
Other important issues.

Help?
Thanks.


texts: Label only showing when polygon is in edit mode

texts:
I am trying to label multiple parts of a polygon.

As this is not working, I decide to select the part not being labelled, make a new polygon of that ('sss'; magenta) and label that single part polygon (intending to make the layer invisible).

However the magenta polygon's label '18' only appears when the polygon is in edit mode and cannot be seen when the polygon is saved (although the polygon itself is still visible).

Also, I do not have 'Only draw labels which fit completely within feature' ticked.
And do have 'Show all labels for this layer' ticked.

Does any one why a label might only be visible in edit mode?


texts: Displaying polygons with colour gradients that changes at value 0 using QGIS?

texts:
I am plotting some model outputs in QGIS and let's say I have got a given number of polygons and for each one a value which can be negative OR positive.

I need to assign to each polygon's value a colour and it needs to change between positive and negative, e.g. red for negative and green for positive.

I think that the Graduated option is the one that best suits my needs but I don't know how?

Also if this can be done, how can I apply it automatically to other layers?  I have quite a few maps to do, and to do it manually would take a long time.

The idea is something like this but with polygons:


texts: Where can I go to have a rant about Open GEO Data?

texts:
Which is the most popular GIS forum or website that's not under the control of our corporate or government overlords?

This place is great for questions and answers, but if I want to have a discussion or in this case a rant about the state of open data access in the GIS community, which is the most popular forum for such things?

p.s Don't worry, I'll be marking this as a community wiki myself ;-)


texts: What I know is that WMS is a tiled image of the layer, so my question is how does it get the attributes of the layer for queries and info-windows as long as it's an image? and why can't we edit them then?


texts: Updating file locations in QGIS projects?

texts:
I have about 30 different QGIS Project files for 30 different clients. Each is saved onto my computer and all their relevant vector and raster files are all saved locally.

We are now in the process of moving office and upgrading our computers.

As a result, I need to move all my files from my computer onto the network drive, thus destroying any previous file references.

Aside from manually clicking on each file and looking for each ones new location, is there any other way around this.

(Each project could have a few hundred links to update)


texts: Adding unique id automatically by field calculator as counter in QGIS 3

texts:
I tried a lot by field calculator but it gives me null as you show in table and another try it gives me value 1 in all rows.

Most of examples, on internet they are exlaining old version and I have QGIS version 3  and I am not sure how can I add id for all rows .

If I want to generate random id, how can I do it also.

i did as suggested in answer  $id

but could i ask why the values are not  looks like counter ,it looks random .
if i need to increase by one ,how can i do that


texts:
I'm trying to get all the points from a Self-Intersection LineString.

On LineString segments intersection with no vertex on intersection:
ST_Intersection make the job;
Now, how catch the vertex intersection like image? Just crosses points...


texts: Setting number of Decimals for QgsPalLayerSettings object

texts:
How can I set the number of decimals for QgsPalLayerSettings object? I used the following code (from QGIS Python Programming Cookbook).e

label = QgsPalLayerSettings()
label.readFromLayer(shp)
label.enabled = True
label.fieldName='tmp_Averag'
label.placement= QgsPalLayerSettings.AroundPoint
label.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'8','')
#label.decimals = 2
label.writeToLayer(shp)
QgsMapLayerRegistry.instance().addMapLayers([shp])

I tried to read the API document and I thought

label.decimals = 2

does the job. But it did not work.


texts: Legend and coloring of markers with plotGoogleMaps in R?

texts:
I have a data frame (X) that looks like this:

code   cluster    name   longt      latit
101      1         A   -89.6171 35.24992
105      1         B   -90.0154 35.10510
106      2         C   -89.5213 34.93277
111      1         D   -86.8642 36.34807
113      3         E   -86.6125 36.19958
114      2         F   -90.4611 30.02196

I would like to plot the observations based on the cluster using plotGoogleMaps function. I have done the following:

X1 <- X
coordinates(X1) <- ~longt+latit # convert to SPDF
proj4string(X1) <-  CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
ic <- iconlabels(X1$cluster, height=10,colPalette=rainbow(4), icon=TRUE)
m <- plotGoogleMaps(X1,filename='myMap.htm', iconMarker=ic,
    colPalette=rainbow(4),mapTypeId='ROADMAP',layerName = 'Clusters')

But I am facing two problems:

when I use rainbow(3) as colors, I get only 2 colors in my legend (instead of 3). But when I use rainbow(4), I get 4 colors.
The legend doesn't represent the clusters, i.e. I don't get a color for each cluster (like one color for #1, one color for #2, one color for #3,...)
The Icons plotted on google map are only 2 colors (Green and Red). Clusters 1 and 2 are all red, and cluster 3 is green. I would like to obtain a distinct color for each cluster.


texts: Measure of isolation (better than distance to the closest)

texts:
I have square locations (see image). For my SDM, I need some measure of isolation of a square (related to the other squares or a particular specified set of squares). The most simple one would be to take the distance to the closest square, but that's not the best, because it might be like 20 other squares at that distance, or it might be just that single square and then vast area of nothing around... Taking just like average distance of 5 closest squares seems kind of hoky poky... Isn't there a better and statistically more robust measure of isolation?

The purpose of this measure: I want to see the distribution of isolation of the squares which need to be predicted with the SDM and compare it to the distribution of the isolation of the clustered cross-validation folds (in the picture), to see if the characteristics of isolation are similar.


texts:
I am searching for a gis/map engine SDK for WPF, that supports loading and editing CAD drawings onto a map (edit = switch layers on/off), with an option to use only local/network resources. Our users are diverse and some might use an actual map (eg bing maps), others might only need CAD drawings. The GIS part seems to be less important right now, but we want a future proof solution, and GIS will definately become a high priority over time.

So far I have been looking mostly at ArcGIS, because it seems to be able to do a lot. But then I discover I am using an SDK that is about to be retired, and the new one does not support WPF controls on the map itself. No good. The map symbols need to very interactive and configurable at runtime. Maybe there is another way to achive this in ArcGIS? Also, neither of the SDKs seem to be able to load CAD drawings, unless they have already been converted to ArcGIS layers in ArcMap.

Other candidates include MapDotNet and ThinkGeo, but nothing really seems optimal so far.

Edit: If anybody is curious, we ended up choosing ThinkGeo.


texts: How can I assign a Vehicle Routing Problem to visit orders multiple times?

texts:
I am developing a new set of routes for my organization facility inspectors to minimize travel time. I have used ArcMap's Network Analyst to solve this for seven routes (one for each day of the week), and it worked well. However, after presenting it to the end users, I learned that different facilities are on different inspection cycles- some are visited weekly, some daily, and some two, three, or four times per week.

Is there a way I can have my orders be revisited?

I've thought of copying those orders that need multiple visits and then assigning them to separate routes or providing them with different time-windows, but that would make the assignments arbitrary and defeat the point of the analysis.

I'm running ArcMap 10.2.2


texts:
I am using following code to display geotiff file. Its only a small file 1,5MB but when I am moving map I have online white screen which gets loaded after I finish moving but it is a bit slow for me. I am scared how slow will it be with 1,2GB map. So is there a way to make this faster or to render map just as it is showed (not wait until user finishes moving)

Here is the code:

public class App
{
    private StyleFactory sf = CommonFactoryFinder.getStyleFactory();
    private FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
    private GridCoverage2DReader reader;

    private JMapFrame frame;
    public static void main( String[] args )
    {
        App me = new App();
        try {
            me.showMap();
        }catch(Exception e){
            e.printStackTrace();
        }
    }

    public void showMap() throws IOException {
        File rasterFile = new File("C:/Users/Daniel/Desktop/release-1800-x64-gdal-1-11-1-mapserver-6-4-1/out.tif");
        AbstractGridFormat format = GridFormatFinder.findFormat(rasterFile);
        reader = format.getReader(rasterFile);

        Style rasterStyle = createRGBStyle();
        Layer rasterLayer = new GridReaderLayer(reader, rasterStyle);

        final MapContent map = new MapContent();
        map.setTitle("ImageLab");
        map.addLayer(rasterLayer);

        frame = new JMapFrame(map);
        frame.setSize(800, 600);
        frame.enableStatusBar(true);

        frame.enableToolBar(true);

        frame.setVisible(true);
    }
    public Style createRGBStyle() {
        GridCoverage2D cov = null;
        try {
            cov = reader.read(null);
        } catch (IOException giveUp) {
            throw new RuntimeException(giveUp);
        }
        // We need at least three bands to create an RGB style
        int numBands = cov.getNumSampleDimensions();
        if (numBands < 3) {
            return null;
        }
        // Get the names of the bands
        String[] sampleDimensionNames = new String[numBands];
        for (int i = 0; i < numBands; i++) {
            GridSampleDimension dim = cov.getSampleDimension(i);
            sampleDimensionNames[i] = dim.getDescription().toString();
        }
        final int RED = 0, GREEN = 1, BLUE = 2;
        int[] channelNum = { -1, -1, -1 };
        // We examine the band names looking for "red...", "green...", "blue...".
        // Note that the channel numbers we record are indexed from 1, not 0.
        for (int i = 0; i < numBands; i++) {
            String name = sampleDimensionNames[i].toLowerCase();
            if (name != null) {
                if (name.matches("red.*")) {
                    channelNum[RED] = i + 1;
                } else if (name.matches("green.*")) {
                    channelNum[GREEN] = i + 1;
                } else if (name.matches("blue.*")) {
                    channelNum[BLUE] = i + 1;
                }
            }
        }
        // If we didn't find named bands "red...", "green...", "blue..."
        // we fall back to using the first three bands in order
        if (channelNum[RED] < 0 || channelNum[GREEN] < 0 || channelNum[BLUE] < 0) {
            channelNum[RED] = 1;
            channelNum[GREEN] = 2;
            channelNum[BLUE] = 3;
        }
        // Now we create a RasterSymbolizer using the selected channels
        SelectedChannelType[] sct = new SelectedChannelType[cov.getNumSampleDimensions()];
        ContrastEnhancement ce = sf.contrastEnhancement(ff.literal(1.0), ContrastMethod.NORMALIZE);
        for (int i = 0; i < 3; i++) {
            sct[i] = sf.createSelectedChannelType(String.valueOf(channelNum[i]), ce);
        }
        RasterSymbolizer sym = sf.getDefaultRasterSymbolizer();
        ChannelSelection sel = sf.channelSelection(sct[RED], sct[GREEN], sct[BLUE]);
        sym.setChannelSelection(sel);

        return SLD.wrapSymbolizers(sym);
    }

}


texts:
I need to learn about the WFS Filter.

is their any website or book where we can learn about WFS and its filter by example.

i have seen the OGC website where we have only xsd.


texts:
I have a question related to the Update method, in case the feature doesn't exists, the update method will generate a new feature? (i.e. will behave as Insert)?

Or maybe, in contrast, an insertion using an existing fid (feature id) will perform an update over this feature?

Thanks in advance.


texts: How to add other column into pgrouting query result in pgAdmin output pane?

texts: I'm working on my final research project now and I'm a new user in pgrouting. Is there any chances to add another column name to display in the output pane? I'd like to add road name column so that I can easily know the road, for the final result I'll build a web app


texts: Exporting attribute table with 100,000+ features into Excel file using ArcGIS Explorer?

texts:
I need to export the underlying data contained in these GIS files (http://wildfire.cr.usgs.gov/firehistory/data.html) into a .txt or .ascii file.

I am using ArcGIS explorer so I do not have the "Table to Excel function."

Some of the layers have less than 100,000 features so I've been creating attribute tables and then exporting those as .txt files. Some of the layers have >100,000 features. I'm at a loss for how to export those into an excel file.


texts: Can ArcGIS Server WMS display a graphicslayer?

texts:
I have tried to display a graphics layer with the WMS functionality of ArcGIS Server, but I cannot see the graphics in the image of the getmap response.

I have done is manually by publishing a service with a graphics on an annotation layer, associated this annotation layer with a specific feature layer, and published it, but the graphics is not appearing, only the features.

And I have done it programmatically by creating a graphicslayer and assign it to a layer, and its the same, graphics are not appearing on the map image returned by WMS.

When I do a REST call to the service I can see the graphics in both cases.

Are graphics not supported with WMS or am I doing something wrong?

I am using ArcGIS Server 9.3.1, Java version.

Nicole


texts: Checking if feature visible using PyQGIS

texts:
I am developing a QGIS plugin in python.

I need to allow the user to select multiple features of the current active layer using an extent.

I am using this code to do that:

self.iface.activeLayer().selectByRect(extent)
selectedFeaturesList = list(self.iface.activeLayer().selectedFeatures())

The problem is my layer is categorized depending on an attribute, and sometimes the user shows/hides some of these categories. And I want the plugin to select only visible features.

So what I want is something like:

for feat in selectedFeaturesList:
    if feat.isVisible():
        visibleSelectedFeatures.append(feat)

I have searched a lot and couldn't find any way to determine the QgsFeature visibility.


texts: QGIS Crashes on USB/Storage device Plug-in/out

texts:
I am using QGIS Version 3.4.0-madeira on Windows 10, and my QGIS crashes each time, a USB device is plugged in/plugged out, or even any network drive gets disconnected with while running QGIS. This only happens when QGIS is open, and a storage drive is connected.

It also makes my QGIS Slow, when the network is busy while copying to network drive or USB.

I read it on QGIS Forum, that this is a bug in the QGIS 3.4. How can I fix this issue ?


texts: Transform CRS for all geometry objects in SQL-Database in one step

texts: I would like to easily change the CRS for all my geometry objects in a PostGIS/SpaitaLite/GPKG database in one step. In my case, all objects are in EPSG:4326 and I want them in EPSG:3044


texts:
I am patching a few rasters together. For my first raster I would like to create a buffer around it and exclude data from the next raster when patching. I need to interpolate between the first two rasters, because this area in the second raster is suspicious.

I created a buffer layer like:

r.buffer --overwrite input=raster1 output=raster1_buffers distances=x1,x2,x3,x4

I then use mapcalc to create a new raster which has null values within the buffer and raster2's values outside it.

r.mapcalc "raster3 = if(raster1_buffers<x1, null(), raster2)"

This returns all -nan values. However, if I swap the position of null() and raster2, mapcalc will return a raster with raster2's values within the buffer zone and nan outside.

Can someone shed some light on this problem?


texts:
I'm running QGIS Las Palmas on Mac OS Mavericks but after enabling all features in the Manage Plugins menu, I can't find the plugin I want.
I downloaded it from the internet and wish to install it manually. From what I understood I should be able to reach the plugins folder in qgis/contents/python/plugins but i don't have a python folder when I open 'contents'.
All I have is qgis/contents/plugins and then a bunch of other folders, but no python.

Can anyone guide me through manually installing the plugin on the correct folder, please?


texts: openlayers 3 ol.source.ImageWMS.getGetFeatureInfoUrl returns bad XY

texts:
I have a singleTile layer created using ol.source.ImageWms and ol.layer.Image. Now I am trying to getFeatureInfo on this layer and no matter where I click on the map, the coordinates are X=50, Y=50

I looked into code of getFeatureInfoUrl:

var extent = ol.extent.getForViewAndSize(
                evt.coordinate, viewResolution, 0,
                ol.source.ImageWMS.GETFEATUREINFO_IMAGE_SIZE_);

var x = Math.floor((evt.coordinate[0] - extent[0]) / viewResolution);
var y = Math.floor((extent[3] - evt.coordinate[1]) / viewResolution);

and for my parameters it does always calculate x,y as 50, 50 for default GETFEATUREINFO_IMAGE_SIZE_=[101,101]. With the same parameters ol.source.TileWms.getGetFeatureInfoUrl returns correct coodinates but the calculations there are different.

Is this a bug or I am doing something wrong?


texts:
I want to apply a symbology scheme for the final layer of my model.
Rules are -

All negative - Symbol A
0 - Symbol B
All positive symbol C
unknown - symbol D

I am not sure how I can define these intervals under Classification of Layer Properties. The minimum and maximum can vary since it will be used in a model. I tried to define manually (among other options - defined, standard deviation...etc.) but it worked only for that sample and failed for another example. I am using 'Apply Symbology from Layer' to display the model output finally.


texts: Output of a value at a given date and location in Google Earth Engine

texts:
I calculated the Albedo of a Greenland and would like to code it in a way that GEE gives out the value of the Albedo for a given date and location or to have an overview in a table. So far I was only able to implement a time range. I am unsure on what to code so that I receive the value of the location (e.g. with giving coordinates)
My Albedo calculation code:
var albedo = function(image){
  var alb = image.expression(
  "((0.356*blue)+(0.130*red)+(0.373*nir)+(0.085*swir)+(0.072*swir2)- 0.018)/ 1.016",
  {
    'red': image.select('B3'),
    'blue':  image.select('B1'),
    'nir':  image.select('B4'),
    'swir':  image.select('B5'),
    'swir2':  image.select('B7')
  });
  return(image.addBands(alb.rename("albedo")));
};

var dataset = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA")
                  .filterDate('2020-06-01', '2020-06-30')
                  .filterBounds(Greenland)
                  .map(albedo);

print("dataset",dataset);

var myAlbedo = dataset.select("albedo");
print("myAlbedo",myAlbedo);


texts: Southern Newfoundland and New Brunswick buildings data

texts:
I am looking for GIS data for southern Newfoundland and New Brunswick.

For NFLD: - buildings, property line/polygon, basic contours, topographic items
For NB: - buildings

Let me know of any data sources or websites that may be of interest.


texts: Seeing basemap created by Oracle map builder in database

texts:
I am new in oracle, I tried to create a style, then different themes and then a base map. I saved all my work, but I can't see the base map in the database.

Did I miss something?

The interesting part is that I have a spatial table, which has a geometry field, but I can't see that in Map Builder.


texts:
I want to import a .las file into QGIS 2.0.1. My guess is that I need to convert it into a format that is readable in QGIS. But how?

I tried using LASTools. I downloaded it, but I don't understand how to use it? Any ideas?


texts: QGis cannot retrieve a text value from QLineEdit

texts:
I am a developing a small plugin with python for QGis. I have a QLineEdit widget
and I need to take the input from there. The input is a 8 character number string.
However when I try to retrieve the text from the QLineEdit python returns a unicode object which does not seem to have any string stored in it.

The QLineEdit field is declared in the constructor my form as follows

self.combo = QtGui.QLineEdit()

The text is retireved during a signal-slot call as follows

the_text = self.combo.text()
print the_text
print len(the_text)

This should be very simple. Am I doing something wrong or this is a QGis Bug?


texts: Excluding feature classes from list using ArcPy

texts:
I'm new to Python and trying to buffer some feature classes within a list in a Geodatabase. I used the Clip* wildcard to select the files, but I want to exclude three of those files, so I don't use the buffer on them. I tried the if not fc.startswith('Clip_blabla') but it won't run, as it says invalid syntax. Any ideas?

fcList = [arcpy.ListFeatureClasses("Clip*") if not fc.startswith('Clip_blabla')]

#Create a loop to buffer Clip*
bufferList = []

for fc in fcList:
    arcpy.Buffer_analysis(fc, fc + "Buff", "700 Meters", "", "", "ALL")
    #Add each buffer feature class to a list of feature classes
    bufferList.append(fc + "Buff")
    print(fc + "Buff")
    if fc == "ClipProjMelb_UrbanRes":
        break


texts: I have a bunch of telephone numbers (just the area codes), and I'd like to geolocate them and make a map, but I haven't had too much luck googling for shapefiles for Area Codes. Surely these are in the public domain?


texts: Calculating Date Field with today's date from datetime.date.today()?

texts:
I'm working on part of tool that will use arcpy.CalculateField_management to add the current date to the attribut table. I've wandered far and wide on the interwebs and can't seem to find the resolution to this issue.

When using this code, i get the value "12:00:00 AM"

    input = r'C:\GIS\CARGIS\SHAPES.gdb\CRASH_ON_2013'
    today = datetime.date.today()
    dte = today.strftime('%m/%d/%Y')
    arcpy.CalculateField_management(input,"DTCARXTRCT",dte,"PYTHON")

When using this code, I get the value "06/19/1905" from today's date of "10/07/2014"

    input = r'C:\GIS\CARGIS\SHAPES.gdb\CRASH_ON_2013'
    today = datetime.date.today()
    dte = str(today)
    arcpy.CalculateField_management(input,"DTCARXTRCT",dte,"PYTHON")

Any thoughts on what I'm doing wrong? I'd like to avoid using the update cursor, but will do so if it's the last option.


texts: rendering a single band raster - new maxvalue

texts:
I have a downloaded raster file, singlebandgray, and want to redefine it's maximum value to a higher integer value. Via the QGIS desktop it works smoothly by hand via the Properties.
Via Python I was not able to figure it out. I tried something like this (some snippets that I found somewhere)

fcn = QgsColorRampShader()
fcn.setMaximumValue(66)
shader = QgsRasterShader()shader.setRasterShaderFunction(fcn)
renderer = QgsSingleBandGrayRenderer(rasterlayer.dataProvider(), 1, shader)
rasterlayer.setRenderer(renderer)

The raster file is saved as "rasterlayer". Its minimum value is 0. The maximum value should be 66 instead of a lower integer. Everything else should stay equal.

I should add that I am not looking for a solution which loads the style with an external .qml file. I tried that and it works fine. But it is not my desired solution.


texts: Change the hover contents of the control navigation buttons in Mapbox GL

texts:
On my Mapbox map, I added the navigation control elements proposed by the Mapbox API (zoom and rotation).

map.addControl(new mapboxgl.NavigationControl());

When hovering over these buttons, a content appears as "Zoom in".

This content is declared in the html code of Mapbox, it corresponds to the "title" value.

In the context of a French web application, I'm trying to translate this content. So, I modify the css code in the following way but nothing happens:

.mapboxgl-ctrl-zoom-in:hover {
  content: "test";
}

I don't know if it's possible to change the content directly. I think that an alternative solution would be not to add these navigation elements in the javascript code but directly by declaring them in the html code in order to be able to modify this value.

You can try by clicking on the following link :
JSFiddle Demo


texts:
I'd like to know how to adjust properties of ol/Kinetic. And also, what are the default properties when using the default ol/interaction/DragRotateAndZoom ?

Properties are decay, minVelocity and delay. I tried to get a sense of those values by tweaking them quite randomly but nothing seems to create a panning with a nice inertia. I also tried to find answers inside the source code here, but I couldn't find any.

Here's my attempt

var interactionPan = new ol.interaction.DragPan ({
    kinetic: new ol.Kinetic(-20, 1, 30)
});


texts: I'm trying to display a rather complex shapefile in a browser. I'm using QGIS's QGIS2Web plugin. The shapefile is only about 7 mg, but it is made up of about 66,000 line segments, so it isn't rendering properly. What are my options? Is there some way to only render parts of the geojson at lower zoom levels, so the browser doesn't try to render the whole thing?


texts:
I want to know the difference about TOA T1 and T2. Are they two repeating collections?  I want to use the imagecollection (merge from Landsat TOA T1 and T2).
How can i achieve it?


texts: Updating GeoTools 18.1 to 22.4 migration issue

texts:
My working code using GeoTools 18.1 is below:
for example:
com.vividsolutions.jts.geom.Geometry out = pGeom.buffer(distance.doubleValue(unit));
com.vividsolutions.jts.geom.Geometry retGeom = out;
// reproject the geometry to the original projection
if (!(origCRS instanceof ProjectedCRS)) {

    retGeom = JTS.transform(out, fromTransform);

  } catch (MismatchedDimensionException e) {
    System.err.println(e.getMessage());
  } catch (org.opengis.referencing.operation.TransformException e) {
      System.err.println(e.getMessage());
}

The issue is JTS.transform(...), expects something called an Envelope and not a Geometry.  How do I convert a Gemoetry into an Envelope class for the latest GeoTools?


texts:
My question is related to
Why is Shapely not installing correctly?

However it is not working and I am struggling to find information/understand how to install it.
I am running my own environment using venv.

I have downloaded the binary correct to my Windows architecture (64-bit) and python version (3.7.3) from https://www.lfd.uci.edu/~gohlke/pythonlibs/#shapely.
My pip version is updated to last version (pip-19.0.3).
I have placed the file in the root of my app where I am running my venv from.

In my cmd terminal I am trying

pip install Shapely-1.6.4.post2-cp37-cp37m-win_amd64.whl

But this just give me

ERROR: Shapely-1.6.4.post2-cp37-cp37m-win_amd64.whl is not a supported wheel on this platform.

What am I doing wrong?


texts:
This may well be too far off topic, but I'm sure some here will use the QGIS Redmine issue tracking system. I've recently come to register an issue but keep getting the following issue.

Am I being daft or is there really no box marked 'Category' on the screen? In fact, there's no box left blank at all except assignee.

Does anyone use the issue tracking system, and if so can you let me know what I'm doing wrong here?


texts:
I have used Make XY Event Layer  to Creates a new point feature layer based on x- and y-coordinates using python and it worked for me. Now as a need of my project i want to apply the same in c# using ARCObject SDK10.0 . Additionally i also want to create a table in geodatabase to hold all the xy points in the table with following field : ObjectID, Xcoor,ycoor, Shape.
I would appreciate any kind of hint how to do this in c#.


texts: Difficulty concerting coordinates system of data frame to raster coordinate system (trasformation) in R

texts:
Can someone help me understand why I have issue getting my dataframe (point coordinates into the same coordinate with my Raster data?

My dataframe with point files (lets call it trainData_crs)

trainData_crs<- read.csv(file.csv)

class(trainData_crs)
[1] "data.frame"

The data frame had no coordinate system

crs(trainData_crs)
[1] NA

summary(trainData_crs)
 PlotType      Class                      Latitude       Longitude
 closedforest:22   Min.   :1.000   Min.   :17.81   Min.   :24.72
 grass       :34   1st Qu.:2.000   1st Qu.:17.86   1st Qu.:24.98
 lowforest   :24   Median :3.000   Median :17.92   Median :25.02
 mediumforest:26   Mean   :2.904   Mean   :17.95   Mean   :25.02
 shrubs      : 4   3rd Qu.:4.000   3rd Qu.:18.01   3rd Qu.:25.10
 unknown     : 3   Max.   :7.000   Max.   :18.46   Max.   :25.16
 water       : 1

Create a TIFF raster and convert it to Rasterbricks

all_landsat_bands<-list.files ("file.tif")

img <- brick(all_landsat_bands)

img <- brick(LandsatAprJune_2019_st)
 class(img)
[1] "RasterBrick"
attr(,"package")
[1] "raster"
 names(img)
[1] "B1"  "B2"  "B3"  "B4"  "B5"  "B6"  "B7"  "B10" "B11"

First convert data frame to Latitude and Longitude coordinate numbers

coordinates(trainData_crs)<-~Latitude +Longitude

proj4string(trainData_crs) = CRS("+init=epsg:4326")

Convert to same projection

trainData_crs_sp <- spTransform(trainData_crs, crs(img))

Confirm if points and image are in same coordinates

crs(trainData_crs_sp)
CRS arguments:
 +proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m
+no_defs

 crs(img)
CRS arguments:
 +proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m
+no_defs

Plot them to  check if the points overlap with rasterBrick

plotRGB(img,
    r = 4, g = 3, b = 2,
    stretch = "lin",
    scale=1000,
    axes = TRUE,
    main = "RGB composite image\n Landsat Bands 7, 3, 2"); plot(trainData_crs_sp, add=T)

But the difficulty that I have is that, even after converting the dataframe of points to rasterbrick coordinates, they still do not fall in the same area or overlap.


texts: Arcpy Update Cursor with Character Length

texts:
I am attempting to create an update cursor to concatenate fields of a parcel ID based on the amount of characters in the MAPNO field.  This is important because the spacing between MAPNO and LOT of the parcel ID is dependent on the number of characters in MAPNO.  Below is my current syntax.  The delete and append work just fine.

I am having an issue with the if len(row[1]) == '2' and the proper syntax to concatenate the Parcel_ID field.

import arcpy

def main():

inputParcels = arcpy.GetParameterAsText(0)
targetParcels = arcpy.GetParameterAsText(1)
parcelFields = ["Prefix","MAPNO","LOT","SUBLOT","Parcel_ID"]

try:
    arcpy.DeleteRows_management(targetParcels)
    arcpy.Append_management(inputParcels,targetParcels,"NO_TEST")

    with arcpy.da.UpdateCursor(targetParcels, parcelFields) as cursor:
        for row in cursor:
            if len(row[1]) == '2':
                row[4] = ((row[0][-2:]) + " " + row[1] + "    " + row[2] + "" + row[3])
            if len(row[1]) == '3':
                row[4] = ((row[0][-2:]) + " " + row[1] + "   " + row[2] + "" + row[3])
            if len(row[1]) == '4':
                row[4] = ((row[0][-2:]) + " " + row[1] + "  " + row[2] + "" + row[3])
            if len(row[1]) == '5':
                row[4] = ((row[0][-2:]) + " " + row[1] + " " + row[2] + "" + row[3])
            cursor.updateRow(row)

    del cursor

except Exception, e:
    import traceback, sys, os
    tb = sys.exc_info()[2]
    arcpy.AddError("%s: Line %i"%("filename",tb.tb_lineno))
    arcpy.AddError(e)
finally:
    arcpy.AddMessage("Script is complete")
if __name__ == "__main__":
    main()


texts: Writing a loop in Google Earth Engine and save all files as unique

texts:
https://developers.google.com/earth-engine/datasets/catalog/NASA_NEX-DCP30_ENSEMBLE_STATS

From the above link I need to download daily data for two variables (pr, tasmax) for two time periods: 1980-2005 and 2040-2060 for two sites: site a (lon = 83.1499, lat = 23.14985) and site b (79.477034, 11.015846)

Using the link below I managed to do it for a single site X variable X date combination as shown

Download NASA Earth Exchange Global Daily Downscaled Projections for single lat lon

var startDate = ee.Date('1980-01-01');
var endDate = ee.Date('2006-01-01');

// get the dataset between date range and extract band on interest
var dataset = ee.ImageCollection('NASA/NEX-GDDP')
              .filter(ee.Filter.date(startDate,endDate));
var maximumAirTemperature = dataset.select('tasmax');

// get projection information
var proj = maximumAirTemperature.first().projection();

var point = ee.Geometry.Point([83.1499, 23.14985]);

// calculate number of days to map and extract data for
var n = endDate.difference(startDate,'day').subtract(1);

// map over each date and extract all climate model values
var timeseries = ee.FeatureCollection(
ee.List.sequence(0,n).map(function(i){
var t1 = startDate.advance(i,'day');
var t2 = t1.advance(1,'day');
var feature = ee.Feature(point);
var dailyColl = maximumAirTemperature.filterDate(t1, t2);
var dailyImg = dailyColl.toBands();
// rename bands to handle different names by date
var bands = dailyImg.bandNames();
var renamed = bands.map(function(b){
var split = ee.String(b).split('_');
return ee.String(split.get(0)).cat('_').cat(ee.String(split.get(1)));
 });
// extract the data for the day and add time information
var dict = dailyImg.rename(renamed).reduceRegion({
reducer: ee.Reducer.mean(),
  geometry: point,
  scale: proj.nominalScale()
}).combine(
  ee.Dictionary({'system:time_start':t1.millis(),'isodate':t1.format('YYYY-MM-dd')})
);
return ee.Feature(point,dict);
})
);

Map.addLayer(point);
Map.centerObject(point,6);

// export feature collection to CSV
Export.table.toDrive({
  collection: timeseries,
  description: 'a_hist_tmax',
  fileFormat: 'CSV',
 });

How do I create a loop to process all of the combinations (each site X variable X dateset combination) and save them as unique files in my drive


texts: Choosing URL to use arcpy.SignIntoPortal_server() correctly?

texts:
I'm trying to write some code to automate updating a feature service. I found some code online that uses,

arcpy.SignIntoPortal_service(username,password,portal)

to sign into the portal service.

I looked up the documentation of the code here. The documentation states that the portal parameter doesn't actually do anything and that,

  "[the] default value is URL for the ArcGIS portal currently chosen by
  the user in Desktop Administrator."

Should I add the rest services url to the list in the desktop administrator or should I add the admin url?


texts: Adjusting edges around markers in legends in GeoPandas plot

texts:
I am trying to adjust the edges around the markers in my legend.
Currently my code looks like this

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(15,10))
lri_map.plot(column= '30',
           ax= ax,
            alpha=1,
           scheme = 'NaturalBreaks',
             cmap = 'Reds',
             edgecolor = 'black',
            linewidth = 0.5,

           missing_kwds = {
               'color' : 'lightgrey',
               'hatch': '////' ,
               'label' : 'Missing Values'

           },

           legend=True,
           legend_kwds={'loc': (0.05, 0.2), 'fontsize' : 10, 'facecolor' : 'white', 'edgecolor': 'white' },

           )

ax.set_axis_off()

Where lrimap is a GeoDataFrame.
I produce the following graph:

I would love to have an edge (line) around the markers in the legend, however it seems that by default this is turned off. Is there any way to add a small black line around the markers (it becomes otherwise hard to see the light red marker)? Somehow it does work for the missing data marker...


texts:
I attempted to use GDAL without installing due to not having rights to install on this computer. I managed to get several gdal modules working, but when I attempted to use ReadAsArray() I received the import error

ImportError: No module named _gdal_array

I looked through the python folder and found that gdal and gdalconst both had pyd files and compiled files, but gdal_array does not.

When I ran VersionInfo() I got the value 1910 for gdal if that helps. I also have python 2.6 and numpy

Does anyone know a way I could fix this problem without using a gdal installer?


texts:
My following query gives output of 72 rows which is correct.

SELECT distinct on (st.gid) sv.timing, st.gid, st.area
FROM servers sv
INNER JOIN sites st
ON (st_distance(sv.geom, st.geom)<8.5) where
        sv.workmode = 1 AND sv.timing >= timestamp '2018-02-02 00:00:00'
        and sv.timing <= timestamp '2018-02-02 23:59:59';

But if I remove the condition on timing and want to find out whole output, it misses one row which is between this time. Why is it happening? What is wrong with my query?

SELECT distinct on (st.gid) sv.timing, st.gid, st.area
    FROM servers sv
    INNER JOIN sites st
    ON (st_distance(sv.geom, st.geom)<8.5) where
            sv.workmode = 1;


texts: 48 bit GeoTIFF MapServer 7 WMS not working

texts:
I am having some trouble trying to serve 48 bit GeoTIFF images with a MapServer 7 WMS.

I created two layers for speed purposes with scale limiters.

The low resolution (lowres) layer is served well and looks good, it is 24 bit GeoTIFF data (8 bits per channel). But if the layer changes to the high resolution (highres), 48 bit GeoTIFF data (16 bits per channel), the data is wrong.

I searched for some hints regarding a similar issue but could not find something.

This is a screenshot of the loaded highres layer. It does not matter if I load a bbox inside the browser, QGIS or another WMS client. The highres data looks always the same.

image loaded from highres layer:

https://i.imgur.com/PJVnZW6.jpg

The Mapfile:

MAP
NAME TESTMAP
STATUS ON
SIZE 4000 4000
EXTENT 606000 5400000 642000 5878000

OUTPUTFORMAT
 NAME "GEOTIFF"
 DRIVER "GDAL/GTiff"
 MIMETYPE "image/tiff"
 IMAGEMODE RGB
 EXTENSION "tif"
END # outputformat end

WEB
    IMAGEPATH "/maps/2016"
    IMAGEURL "/maps/2016"

    METADATA
      "wms_title" "TEST"
      "wms_onlineresource"  "http://192.168.3.235:8080/?map=eo" ### recommended
      "wms_srs" "epsg:32718"
      "wms_enable_request" "*"
      "wms_formatlist" "image/jpeg, image/geotiff"
      "wms_format" "image/geotiff"

    END #METADATA end
END #web end

PROJECTION
"init=epsg:32718"

END #PROJECTION end
  LAYER
      NAME EO2016_Full
      GROUP EOGroup
      STATUS ON
      TYPE RASTER
      TILEINDEX "/maps/SHP/2016higres.shp"
      TILEITEM "Location"
      MINSCALEDENOM 1
      MAXSCALEDENOM 20000

      METADATA
      "wms_title"  "HighRes"
      "wcs_extent" "606000 5400000 642000 5878000"
      "wms_srs" "epsg:32718"
      "wcs_formats" "GEOTIFF"
      "wcs_bandcount" "3"
      END #METADATA end

      PROJECTION
      "init=epsg:32718"
      END #PROJECTION end
  END #layer END

  LAYER
      NAME EO2016_Overview
      GROUP EOGroup
      STATUS ON
      TYPE RASTER
      TILEINDEX "/maps/SHP/2016lowres.shp"
      TILEITEM "Location"
      MINSCALEDENOM 20000

      METADATA
      "wms_title"  "LowRes"
      "wcs_extent" "606000 5400000 642000 5878000"
      "wms_srs" "epsg:32718"
      "wcs_formats" "GEOTIFF"
      "wcs_bandcount" "3"

      END # METADATA end

      PROJECTION
      "init=epsg:32718"
      END

  END #layer end

END #map end

Geotiff gdalinfo for one file:

 Driver: GTiff/GeoTIFF
 Files: XXXXX.tif
 Size is 4000, 4000
 Coordinate System is:
 PROJCS["WGS 84 / UTM zone 18S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235604902,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-75],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32718"]]
 Origin = (896000.000000000000000,6144000.000000000000000)
 Pixel Size = (0.500000000000000,-0.500000000000000)
 Metadata:
  AREA_OR_POINT=Area
 Image Structure Metadata:
  INTERLEAVE=PIXEL
 Corner Coordinates:
 Upper Left  (  896000.000, 6144000.000) ( 70d40'26.07"W, 34d46' 9.64"S)
 Lower Left  (  896000.000, 6142000.000) ( 70d40'22.69"W, 34d47'14.39"S)
 Upper Right (  898000.000, 6144000.000) ( 70d39' 7.62"W, 34d46' 6.84"S)
 Lower Right (  898000.000, 6142000.000) ( 70d39' 4.21"W, 34d47'11.58"S)
 Center      (  897000.000, 6143000.000) ( 70d39'45.15"W, 34d46'40.62"S)
 Band 1 Block=4000x1 Type=UInt16, ColorInterp=Red
  Overviews: 2000x2000, 1000x1000, 500x500, 250x250, 125x125, 63x63, 32x32,
 16x16, 8x8
 Band 2 Block=4000x1 Type=UInt16, ColorInterp=Green
  Overviews: 2000x2000, 1000x1000, 500x500, 250x250, 125x125, 63x63, 32x32,
 16x16, 8x8
 Band 3 Block=4000x1 Type=UInt16, ColorInterp=Blue
  Overviews: 2000x2000, 1000x1000, 500x500, 250x250, 125x125, 63x63, 32x32,
 16x16, 8x8


texts: PostGIS - query to get the whole network covered, without overlaps

texts:
I have a situation where I get routes from drivers and store them in a PostGIS database. What I want is to get a GeoJSON from the whole network covered, without overlaps of road segments. In other words, I would like a "cleared" network (see image) of the driven roads.

At first, I thought of creating a trigger for checking if there are any overlaps before inserting a new route. However, a new route could have overlaps itself, so this approach is declined.

Another thought was to split the routes into the minimum, smallest line segments and then delete the ones that are equal to others. However, I think that storing all these segments as different records could be an overkill for the database (correct me If I am wrong).

Ideally, I am searching for a query that would return all network without duplicates (overlaps). Any ideas?


texts: Convert area from square meters to degrees globally

texts:
I have a global raster file with 3600x1800 pixels. I don't have any information about the projection used.

The raster contains one band with a quantity with the unit "per square meter". How can I convert this to be "per square degree", prefereably using GDAL or QGis?


texts: How to create gaussian blur in ArcGIS for Desktop Advanced with Spatial Analyst?

texts:
I am trying to create a gaussian blur for a buffer polygon around a border.

I have Spatial Analyst and ArcGIS for Desktop Advanced, but am not sure if it is possible.

Can anyone suggest a way to do it in ArcGIS for Desktop?


texts: Create a band mask based on pixel locations, not in pixel values

texts:
I have this image loaded in QGIS, with two possible values for the pixels: 1 = orange, 2 = blue

For simplification, assume it represents a province of a country.

I want to extract just the part of the image included inside the limits of the province. Regardless the value of the pixels. I guess I should use a mask for that, but I don't want to mask pixels based on their values. I want to mask the pixels based on their location (inside or outside the limits of the province).

The image is not from a real province. So, I don't know the limits (in that case, I could spatially intersect the raster with the known limits).

One approach would be vectorize this image to get the limits of the fake province as a geometry. That would make things easy. But I'm not sure about how to face that task.

So, suggestions are welcome. Many thanks in advance.


texts: Automatically filling contours by different type of hatching?


texts: Unable to delete attributes from QgsVectorLayer in Qgis 1.8

texts:
I have a layer with several hundreds of attributes in Qgis 1.8. I am trying to remove most of them, but in any case, it also fails with just the first one. I am doing

res = my_layer.dataProvider().deleteAttributes([0])

but res gets the value false if the number of attributes exceeds 403 (a quite extrange limit). Is there a limit of attributes to deal with? How can I overcome this issue?

Update: Example of the issue with three copies of the same layer on Qgis Python console. Same number of features (12), different number of attributes (405, 404 and 403 respectively).


texts: ArcMap Append tool will not copy full contents of input feature class

texts:
I have a versioned SDE geodatabase and and running a series of Appends between feature classes and a stand-alone file geodatabase with a schema that mirrors the input gdb.

I have, for example, an SDE feature class that contains 13k records.  When I append the layer to an empty feature class in the file gdb only about 11k records copy over.  I have recalculated the spatial extents on both.  I even set the environment tool in Append to be the input of both datasets.  It simply won't work and won't copy over the same specific records.  I don't have any subsets of the data pre-selected, fyi.

Can anyone suggest why append is failing to copy over the entire contents of the input layer?


texts: Volume Calculation for non-closed geometries in postgis

texts:
I have a problem regarding the volume-calculation in a specific 3D Dataset.
The data consists of 3D multipatches that represent buildings.

One building consists of one or more non- closed multi-polygons that are open from the inside. (E.g. initially polygons forming a building were closed geometries and overlapped – in the data preparation process the inner-walls had been cut away). As a main drawback volume-calculations with use of the typical functions in postgis is not possible (neither for one part of the building nor for the whole building). The dataset I’m talking about can be found under https://shop.swisstopo.admin.ch/en/products/landscape/build3D2 (additional information/ sample data/ LV03)

In ArcGIS I tried to approximate the volumes with a convex/concave hull (minimum-bounding-volume) neither one did lead to useful results since the deviations from actual volumes were too high. 3D-union could be a possible remedy but since the multi-polygons do no more overlap, the function cannot be applied.

Are there any solutions towards a more accurate approximation of the buildings volumes? (any software tools welcome)


texts: OpenLayers Plugin not working in QGIS 2.8.1

texts:
I am using QGIS 2.8.1. I have installed the OpenLayers Plug, however it is not showing up in my list of plugins. I confirmed it's installed in the plugin manager.

Is anyone else having this issue?


texts: Access-Control-Allow-Origin problem when loading json file from Amazon S3 to display in OpenLayers

texts:
This is my code for displaying geojson in OpenLayers, where "geourl" is directed to a json file in Amazon S3.

geojson_layer = new OpenLayers.Layer.Vector("GeoJSON", {
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
    url:geourl,
    format: new OpenLayers.Format.GeoJSON()})
    });

But I received the error like this:

> Failed to load resource: the server responded with a status of 403
> (Forbidden) "geourl" Failed to load resource: No
> 'Access-Control-Allow-Origin' header is present on the requested
> resource. Origin 'http://localhost:8081' is therefore not allowed
> access. "geourl" XMLHttpRequest cannot load  "geourl" No
> 'Access-Control-Allow-Origin' header is present on the requested
> resource. Origin 'http://localhost:8081' is therefore not allowed
> access.

Nothing changed when I added

self.response.headers.add_header("Access-Control-Allow-Origin", "*")

to my webapp RequestHandler Class.


texts: Python, GeoPandas, how to extract to a list by range of a column

texts:
I have GeoPandas dataframe:

d = {'Angle of point': plot_angles_list, 'Latitude': lat_list, 'Longitude': lon_list}
    df = pd.DataFrame(d)
    gdf = gpd.GeoDataFrame(
        df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude))
    gdf.head()

   Angle of point   Latitude  Longitude                   geometry
0       79.349848  32.959723  35.808097  POINT (35.80810 32.95972)
1       79.349848  32.959723  35.809097  POINT (35.80910 32.95972)
2       79.349848  32.959723  35.810097  POINT (35.81010 32.95972)
3       79.349848  32.959723  35.811097  POINT (35.81110 32.95972)
4       79.349848  32.959723  35.812097  POINT (35.81210 32.95972)

I want to extract to a list only the points that have Angle of point value between 50 - 150, with the expection of 88-92, this values I want to append to another list.


texts: How to calculate multiple travel times using gmapsdistance?

texts:
I am trying to use gmapsdistance on R to calculate multiple transit times from Google Maps.

I have a Google API key and can generate transit travel times for one origin and one destination with the following code:

gmapsdistance("51.7788868858975,-1.48196518842103"
, "51.76732413,-1.233620597" ,"transit",
              "googleAPIKey")

However, i have a list of 1000+ origins and destinations in latitude and longitude that i want to find transit times and can only retrieve one at a time. I've tried to adapt the formula provided here to request multiple requests using GGMAP:

google_results <- rbind.fill(apply(subset(TravelExample, select=c("Origin", "Destination")), 1, function(x)
    gmapsdistance(x2, x1, "transit", "googleAPIKey")))

I get the following error:  Error in gmapsdistance(x2, x1, "transit", "googleAPIkey") :
  Google maps is not able to find a route between origin and destination

I would use GGMAPS, but unfortunately you cannot request transit times from Google using the pakagage.

Any ideas how i can request multiple transit times using gmapsdistance?


texts: How to georeference an old picture through transformation values?

texts:
I got a very old picture of a landscape which I have to georeference the right way. For that, I usually go in my QGIS version 2.12.3-Lyon to “Rasta” and then “Georeferencing”. Normally it’s easy, because I get the coordinates from google maps and then try to find the same points in my picture. Since it is a very old picture, it’s pretty hard to find the right points for the coordinates. I also got some transformation values, but I don’t know how to use them in QGIS and if it is even possible. Someone knows how those values could help me?
I added a few of those transformation values as example.

Edit:
I georeferenced as suggested and added points with some random number and then later on I edited them to the right value. Now I have the problem that my picture is not at all where the points are (if you look very closely you can see a small spot on the left /down). And that the dx and dy values don’t fit with the values I got as transformation value. When I start the transformation I also don’t get the picture in the right place and it is like I didn’t georeferenced it at all… Someone knows how to proceed?


texts:
I tried installing Rtree on Python35 on my windows 7 machine using the pip install command but ran into this error:
OSError: could not find or load spatialindex_c.dll

I also tried this and it gave me the same error:
    C:\Anaconda3\Scripts\easy_install.exe Rtree

I looked up this post on StackExchange:
Installing libspatialindex on Windows 8

However, I am not able to find a binary for Python35. Can you please tell me how to go about installing Rtree?


texts: QGIS Processing - Additional Creation Parameter for GdalWarp

texts:
I am having a hard time getting the "Additional Creation Parameters" to work when using GdalWarp through the processing toolbox.

I am trying to set the nodata value to 0 for the output raster.

processing.runalg("gdalogr:warpreproject","InRaster.tif","EPSG:32736","EPSG:32736",10,0,"-dstnodata 0","OutRaster.tif")

Has anyone used this functionality and managed to get it to work?


texts:
I have looked through several other questions and am still having trouble.  I have a set of 103 polygons that span my entire map.  Each of those polygons is divided into 36 equal sized polygons.  Now I need to iterate through each of those 36 squares, and split them into 100 equal squares.  I've tried getting in contact with the person that split the first 103, but all efforts have failed.  Can anyone point me in the right direction?

Fishnet seems to be out of the question as it spans the entire feature class, not the polygons.

I'm guessing I will be writing a python script.  Any pointers would be appreciated!


texts:
VERSION QGIS 2.1x

i realize plugins to do this actions without use toolbar: copy features / paste to layer / merge selected features => define attribute to use in merge.

So, all is okay but i don't find "merge selected features" in PyQt action or others like :

qgis.utils.iface. ???

What is good expression to call "merge selected features"  action / button ?
All of this actions are in advanced digitize toolbar.


texts: Put the start and the end name of a point object into the layer of the line object that connects them?

texts:
I am new to QGIS. I am given two layers with point objects. In their attribute table each of this objects has a code (name) given previously to it by someone else. Then I import the lines from AutoCAD that connect those points. For those lines I have created 2 fields in their attribute table that say 'Start point' and 'End point'. In the field 'Start point' in the line object table, I should include the code(name) from the first point object, and in the field 'End point' I should include the name of the end point object.

Since there are many lines like this, and each of them connects 2 point objects, is there an option that can do that automatically without me having to type each code(name) differently?


texts: How to combine two string replace functions in one expression?

texts:
I would like to change some of my street labels like "Musterstraße" into "Musterstr." or "Muster Straße" into "Muster Str."

I already did that for the first case:

replace("name", 'straße', 'str.')

and for second case:

replace("name", 'Straße', 'Str.')

But how do I combine both?

I tryed different styles like:

replace("name", 'straße', 'str.') OR replace("name", 'straße', 'str.')

replace(replace("name", 'straße', 'str.')("name", 'straße', 'str.'))


texts: Save virtual fields to a layer definition file

texts:
Is there a way to save virtual fields (with their expression) a layer definition file without them becoming static fields?

I have about a dozen virtual fields in many tables with complex expressions, which need to be replicated in about another 20+ projects.  I would really like to find a way that I can replicate the virtual fields into the other projects without needing to recreate and copy/paste the expressions for every single virtual field.


texts: Why are Character Marker Symbol Dialog not showing up in the Symbol Property Editor dialog?

texts:
I just installed ArcGIS 10.2.1 on a brand new machine and for some reason Character Marker Symbols are not loading into the Symbol Property Editor.

Is anyone experiencing the same issue?  Is there a solution for this?


texts: Extracting Values from Multidimensional Raster

texts:
I am currently doing my dissertation on hurricanes and their relationship with Sea Surface Temperature. I am using a NetCDF file with dimensions: latitude, longitude, and time
I need to extract a value of SST for each point along the hurricane track (therefore matching latitude, longitude and time). This will be for 638 hurricanes, approx. 15-20 points in each, for 40 years.
How can I do this?
I am using ArcGIS Pro 2.6.0


texts: Running local Moran's index in all fields automatically?

texts:
I want to use the local Moran's I index in ArcGIS.

I have a table/layer with 50 columns.

Can I execute moran in all fields of this table automatically or do I have to run the tool 50 times?


texts: I am developing a webapp using opentripplanner. I have followed 2 and five minute tutorials at their websites, there is also a opentripplaner-geocoder.war webapp, but i am unable to figure out that how can i use this geocoder.


texts: How to export a CSV with lat lon WKT from a layer which is not in WGS84?


texts: Does Calculate Geometry tool exist in ArcGIS Pro?

texts: Does ArcGIS Pro have a Calculate Geometry tool available for use?


texts: Accepted data types for wps process JTS:difference in OpenLayers?

texts:
I am trying to use JTS:difference (present as wps service on GeoServer) using OpenLayers API.
In OpenLayers wps-client example,  I have replaced intersection with difference and
I get illegal argument exception.
Can someone please help on what is valid for a: and b: for JTS:difference process?

var map, client, intersect, buffer;

function init() {

    map = new OpenLayers.Map('map', {
        allOverlays: true,
        center: [114, 16],
        zoom: 4,
        layers: [new OpenLayers.Layer.Vector()]
    });

    var features = [new OpenLayers.Feature.Vector(OpenLayers.Geometry.fromWKT(
        'LINESTRING(117 22,112 18,118 13, 115 8)'
    ))];
    var geometry = OpenLayers.Geometry.fromWKT(
        'POLYGON((110 20,120 20,120 10,110 10,110 20),(112 17,118 18,118 16,112 15,112 17))'
    );

    map.baseLayer.addFeatures(features);
    map.baseLayer.addFeatures([new OpenLayers.Feature.Vector(geometry)]);

    client = new OpenLayers.WPSClient({
        servers: {
            opengeo: 'http://demo.boundlessgeo.com/geoserver/wps'
        }
    });

    // Create a process and configure it
    intersect = client.getProcess('opengeo', 'JTS:intersection');
    intersect.configure({
        // spatial input can be a feature or a geometry or an array of
        // features or geometries
        inputs: {
            a: features,
            b: geometry
        }
    });

    // Create another process which chains the previous one and execute it
    buffer = client.getProcess('opengeo', 'JTS:buffer');
    buffer.execute({
        inputs: {
            geom: intersect.output(),
            distance: 1
        },
        success: function(outputs) {
            // outputs.result is a feature or an array of features for spatial
            // processes.
            map.baseLayer.addFeatures(outputs.result);
        }
    });


texts:
How can I create a Digital Surface Model in order to collect the height of buildings and what's exactly needed to create a Digital Surface Model?

Any software provide creating it (ENVI or PCI preferred)

I have a stereo pair satellite images and tried ENVI but I didnt know exactly what I need to have and the requirement to get the perfect DSM

Mostly used software for me is Arcgis and am not used to Imagery


texts: How do I delete part of a polygon (or a number of polygons) outside of another larger polygon in QGIS

texts:
I have got a number of polygons within a single larger polygon. Parts of the smaller polygons go out from under the larger one, and I want to delete the parts of the smaller polygons, that go outside the larger one.
How can this be done?


texts: How to eradicate excess lines in a singly recorded polygonal line (road) network? QGIS 2.8

texts:
I'm a beginner in GIS software.
I'm trying to create a road network shapefile within the reserve I work at. To do this, I've driven each and every road inside the property with my GPS app set to record the track. I've exported this track into QGIS, saved it as an ESRI shapefile, and extracted the nodes. My ultimate aim is to use this line network to create polygons within the boundaries of surrounding and intersecting lines to subdivide the reserve into sections for easier management. You can see what I've done so far below:

My problem is the following -- although I did my best during recording, there were a number of single roads that I recorded more than once (i.e I drove some roads, with the GPS app tracking, 2 or more times each when for example I had to follow a loop road). The result is that instead of having line intersections only at the actual junctions as per satellite/ real-time view, I'm left with this chaotic mish-mash of intersections as you can see below:

This makes it impossible for me to draw out clearly defined polygons from these lines, especially as some lines are just short of intersecting as a result of me pausing and unpausing the track. To add to my woes, the track I recorded has no attribute data (i.e order of nodes) that I can use to break the whole network into segments and connect individually using the 'Join Multiple Lines' function. Because of this, trying to use Geoprocessing functionalities such as 'Dissolve' or 'Eliminate Sliver Polygons' is futile
How do I merge or remove these excess lines and make one completely linear, 2D shapefile of the road network I can use to make clear polygons, considering the above? Please note that I am still using QGIS 2.8.2 and for various reasons am not considering upgrading
I hope I've explained the problem adequately and I have a deadline to follow for this task.


texts: GeoServer create thematic from point to area

texts:
I have some point data like first image.

I want create thematic with this data. But I can't create GeoServer style like this.

You see data on map like this:

How can i create style which i see thematic image on this image? :

I use heat map but heat map does not meet my request:


texts: Adding description to layer using ArcMap?

texts:
This may be quite basic but can't figure it out.

I have a problem when adding a description to a layer (right click layer -> properties -> general: description). When I update the description, then remove the shp and add it again, the description isn't updated. How should I do this to have the updated description saved at the layer's properties? I tried several times with shapefiles and layer files, and it is updated in layer-file properties, but not in shapefile's. I would need the description of a shapefile to be updated because I'm sending shapefiles to another person.

I'm using ArcGIS 10.1 for Desktop.


texts: Removing incorrect (non-ground) points from LiDAR data?

texts:
I have ground lidar points, nonground lidar points and tree lidar points.

The problem is some of the ground lidar points are incorrect. For example, the highest elevation values of tree and building are not over 200 but a few of points in ground lidar points are over 300.

so I want to remove the points higher than 200 and replace them with the surrounding values in ground lidar data.

Any help?


texts: Getting highest resolution data for Costa Rica on Google Earth Engine?

texts:
Unless I am just missing it, I am wondering why the DigitalGlobe WorldView data isn't available through Google Earth Engine even though it is on Google Maps?

I want to get as high resolution data for Costa Rica specifically.


texts: Calling a plugin (quickosm) from a python standalone application?

texts:
I am trying to call the QGIS-plugin quickosm from a standalone Python application and I´m not quite sure how to or if this is even possible.

My QGIS version is 3.10 (using OSGeo4W64) with the quickosm-plugin installed.
Running on Windows.

I tried to achieve this by using QGIS Processing, but I failed.

import sys
import os
import qgis
from qgis.gui import *
from qgis.core import *

from PyQt5.QtCore import *

QgsApplication.setPrefixPath("C:\\OSGeo4W64\\apps\\qgis\\", True)
qgs = QgsApplication([], False)
qgs.initQgis()

qg_project = QgsProject.instance()

from qgis.analysis import QgsNativeAlgorithms
import processing
from processing.core.Processing import Processing

Processing.initialize()

qgs.processingRegistry().addProvider(QgsNativeAlgorithms())

# The extent is calculated somewhere else
alg_params = {'EXTENT': my_extent,
              'KEY': '',
              'SERVER': 'http://www.overpass-api.de/api/interpreter',
              'TIMEOUT': 25,
              'VALUE': ''
              }

feedback = QgsProcessingFeedback()

outputs['BuildQueryInsideAnExtent'] = processing.run('quickosm:buildqueryextent', alg_params, feedback)

qgs.exitQgis()

This fails, stating algorithm quickosm:buildqueryextent not found.
When I try to call QGIS native, GDAL or GRASS algorithms, processing works.

I haven't found a solution to this. Maybe I can tell the script somehow where to find the plugin scripts? Is such a call even possible? Or do I have to do it another way all together?


texts: Normalizing election voter data for choropleth?

texts:
I would like to create a choropleth that shows eligible voters per constituency. This data is in raw count.

How do I normalize the data in this case?

My intention is to show the difference in number of eligible voters for each constituency.

One thing I have done is shown how much each constituency deviates from the average. In some cases this is extreme and I have a range of variation between - 99% and 328%, so also having trouble representing this properly in buckets.


texts:
I am using two Vectorlayers. I use one of them to show Circles with a certain radius. This layer has zIndex 10, the other layer is used to show markers and has zindex 1000.
I want my users to be able to modify the circle, so i defined the following modify:
var modify = new Modify({
            source: this.getCircleLayer().getSource(),
            wrapX: true,
            style: new Style({
                stroke: new Stroke({color:'#AA0000', width:2}),
                fill: new Fill({color: '#AA000059'}),
                image: new CircleStyle({radius:3, fill: new Fill({color: '#AA0000'})}),
                zIndex: 1,
            }),
        })

Now I can see both of my layers and the MarkerLayer is shown on top of the CircleLayer, but the point to modify the circle is shown on top of the Markerlayer and I don't understand why?
Is there a way to change this?


texts: Is there an OGC standard for data synchronization?

texts: ESRI provides two solutions that facilitate synchronization.  They are Editor Tracking for Feature Services (http://server.arcgis.com/en/server/10.3/publish-services/windows/editor-tracking-for-feature-services.htm), which tracks who and when an edit was made, and a built-in way of assigning a UUID to a feature.  Are there any OGC standards which address change tracking in a similar (or better) way?


texts: Extracting raster by mask outputs raster with extent beyond the mask layer's in QGIS

texts:
I'm trying to extract a section of a raster file that overlaps with a polygon shapefile in QGIS 2.18.3. When I run extract, I'm left with a large null value area as part of the extent of the output raster.

I have the polygon (.shp):

And the raster (.tif) which had this large null value area since the first time I clipped one of two DEMs that were merged together to create this layer:

I used Raster -> Extraction -> Clipper tool with the .tif input and .shp mask layer shown above, and it produced this output:

Using clipping option "Crop the extent of the target dataset to the extent of the cutline" didn't do much better:

Here are the settings I'm using in QGIS' clipping tool:

How do I get the DEM clipped exactly to the boundaries of the shapefile of interest in QGIS? Open to recommendations of other open source GIS tools as well but would prefer using this since I'm already up and running and learning my way around it.


texts:
I use ArcGIS 10.1. I have a shapefile with sites and another with the rivers of the area.

I would like to calculate the density of the sites related to the river vicinity. Is it possible? Which type of tool could I use?


texts:
I have done between group pca in R. I want to make an interpolated map using the eigen values of the first three pc axes.

I want to produce a map using each of the three pc axes and finally combine the three maps to produce a final interpolated map.


texts:
I am trying to use HERE application for isochrones calculation. Few months ago everything works fine. I got my APP ID and API key for REST requests and used it in plugin menu. But now something is going wrong. I'm getting Python error:

File "C:/Users/Александр/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\Hqgis\hqgis.py", line 951, in getIsochronesSingle
  ranges = [int(x) for x in intervalArray]
  File "C:/Users/Александр/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\Hqgis\hqgis.py", line 951, in
    ranges = [int(x) for x in intervalArray]
   ValueError: invalid literal for int() with base 10: ''

How can I cope with this error?
I use QGIS 3.4, Win 7-64.


texts: Digitized Shapefile manually using wrong projection: how to fix it?

texts:
I have a raster in correct projection. I used shapefiles, and marked objects (manually), but I found out that my vector layers have different projection, so they have different size and position. Is there a way to correct it using QGIS functional?

That is what I see:

Then, I tried to convert vector projection to raster projection and got that:


texts: Why is my legend in QGIS using cached data?

texts: I have a legend... multiple legends as a matter of fact. Even though I clear them out and start from scratch, I am getting old cached data of what my color ranges used to be, instead of what they've been changed to. I can't get my new colors in the legend... and even the labels are wrong from what I originally had them as.


texts: Spatial and DateTime join in SQL BigQuery

texts:
I have 2 tables in BQ; one of them includes a geography column PointWKT (gps POINT) and second table includes a geography column BoundaryWKT (POLYGONS & MULTIPOLYGONS). Both tables also include DateTime column.
I want to do a spatial join so each PointWKT will be assigned to corresponding BoundaryWKT.
Until today, i was using GeoPandas but now i have to switch to SQL. Is there any way in SQL BigQuery that would substitute geopandas .sjoin like:
join_df = gpd.sjoin(table1, table2, how="inner", op='within')

ANSWER to this question:
SELECT * FROM table1, table2 WHERE ST_Contains (table1.geom, table2.geom)

Another question is, if it is possible to include another condition for joining? Besides the spatial join from geography columns i also need to include join on DateTime column from both tables. Is this even possible? Thanks a lot!


texts: CRS transformation WGS 84 UTM 30N (EPSG: 32630) to ETRS 89 UTM 30N ( EPSG: 25830)

texts:
I have a coordinate table with EPSG 32630 (WGS 84 UTM30N) and when I try to transform them to EPSG 25830 (ETRS 89 UTM30N) the difference between them are marginal.
I´m using QGIS doing the following process: Adding the coordinates as points, then saving them with EPSG 25830 CRS.
I disable the ‘on the fly’ CRS transformation option.
I´m doing something bad? Or is truly marginal this transformation.
Thanks


texts: How can you find ArcGIS version programatically?

texts: Is there a way using ArcObjects.net to find out what version of ArcGIS is installed on a machine (i.e. 9.3., 10.0, 10.1)?


texts: Finding N,S,E,W lat/lon at a determined radial distance from an input lat/long

texts:
I am working on a project where I geocode a users address to get a lat/lon value. With that value, is there a way to find the latitude and longitude of points that are a set distance (eg. 100m) away in a circle?

This diagram might help explain what I am looking for:


texts: How many buildings are there in OpenStreetMap?

texts:
I'm interested in a ballpark estimate of how many buildings there currently are in total in OpenStreetMap.  I'd also like to know roughly how many buildings in total there are in the world, so I can calculate the coverage, but it's mostly the total count that I'm after.

I've attempted to download the Planet.osm data file but the uncompressed version is too large for my computer, so I haven't been able to spin it up and run the analysis myself.

I'm not using any GIS software, but could do so if there's a path to establishing the answer. I was anticipating that the answer could be found using Clojure, Python or by loading the data into a PostGIS instance.


texts: Using Python Parser in ArcGIS Pro Field Calculator?

texts:
I am having a problem in migrating from Vbscript to Python for ArcGIS Pro.

  Building 5, Street 6

What is wrong in this expression?

"{}:{}:{}:{}:{}:{}".format("Building ",!Building_Number!," , ",!Street_Type!), " " ,!Street_Number!


texts: How do I modify the "Playlist" Storymap to include a basemap switching widget?

texts:
I'm fairly new to HTML and the ArcGIS Javascript API. I would like to modify the default "playlist" storymap template to include a basemap gallery widget. Can someone help break down the steps to do this?

This is what the default template looks like:
http://storymaps.esri.com/stories/2013/20towns/

Here is an example of the modified playlist template with the basemap gallery widget.
http://fishandgame.idaho.gov/ifwis/maps/wma/


texts: Transfering files from QGIS to ArcGIS, wrong topology and fantoms objects

texts:
I am working in a project with both QGis and ArcGIS. I do mostly cosmetic work and digitalize with QGIS then transfer the shapefiles on ArcGIS. At work I am the only one to work on QGIS, the other members have ArcGIS, (not our choice, question of tight budget). This also mean I have no support, no one is familiar with the program.

This worked fine as long as I was only transfering ready files to ArcGIS, but lately I worked for a short time on the files in ArcGIS (I projected the files and edited) and eventually I had to take back the files from ArcGIS and edit further with QGIS.

When I send back those files in ArcGIS they are messed up.

The edited shapefiles are right on QGIS (topology fine, TOC fine), but on ArcGIS I have "fantoms" polygons or objects (without any geometry) AND overlapping polygons: the old version of the edited polygon/object plus the new polygon/object.

Though I have a column with unique IDs to help sort that out, some of the fantoms/old polygons get a "0" for ID some other the same ID as the original polygon.

So I delete the mistakes in ArcGIS, the file is fine according to ArcGIS. I double check on QGIS, and surprise, I am now missing polygons there!

This is becoming a big mess, and I am unable to understand where are my mistakes coming from so I can avoid them in the future.

I also have problems with special signs (Swedish signs, ö, ä, å) translating in other signs or letters from on program to the other. I have made some reference excel table to join to the gibberish TOC, but when constantly editing, it is complicated to keep track of the newest version. Any other solution is welcome!


texts: Calculating distances between closest points of each polygon

texts:
I would like to calculate the distances between the closest points of each polygon. I've tried using centroids and both distance matrix and distance to nearest hub, but since I'm relying on centroids, this is associated with some error for my project.

So is it possible to tell QGIS to determine the closest point between each polygon and then calculate it's distance?

QGIS Version: 2.18.18


texts: gdal_translate cannot convert MSG3 channels to TIFF

texts:
following this guide http://www.gdal.org/frmt_msg.html I want to convert my MSG3 HRIT format files in TIFF but I get this error

gdal_translate H-00 0-MSG3__-MSG3________-IR_097___-000008___-201410281200-C_ gg.tif
ERROR 4: `H-000-MSG3__-MSG3________-IR_097___-000008___-201410281200-C_' not recognised as a supported file format.
GDALOpen failed - 4
`H-000-MSG3__-MSG3________-IR_097___-000008___-201410281200-C_' not recognised as a supported file format.

The same error I have on windows.

Can you help me?


texts: How to georeference TIFFs with coordinates from txt file using ArcGIS for Desktop?

texts:
I have a folder with TIFFs and an accompanying .txt file containing coordinates in this format (sample):

    # IMGE       LN   ID             Time         Easting      Northing       Height       Phi      Omega      Kappa
    \00381.tif   2111 21110038       385850.277   32356528.276 5642849.274    1518.200     47.7391   -0.0656    100.1021
    \00391.tif   2111 21110039       385854.176   32356526.170 5643123.992    1517.967     47.7847   -0.1409    100.1607

The images themselves contain no info - the image name is the only connection I have. The images also overlap, so I can't put them together manually, at least not without a lot of work.

Any ideas?


texts: Identifying OpenLayers 3 version at runtime?

texts:
I'm looking for a way to check/confirm which build of OpenLayers 3 is loaded in a browser.  I was expecting that this might be accessible through something like ol.VERSION or ol.getVersion() or similar, but I haven't found anything.

Is there a way to get this information?


texts: ArcPy - Import a directory as a hierarchy of group layers in ArcMap

texts:
I would like to automate the importation of a directory as a hierarchy of group layers in the Table Of Content of ArcMap with Arcpy.
I already saw and adapt these posts and scripts :

Importing file geodatabase (GDB) to table of contents (TOC) of ArcMap?
Adding group layer, renaming it from folder name, then adding layers to it from that folder, using ArcPy/Python?

Unfortunatelly, each code had only one level depth, so there’s no group layer(s) into another group layer. In other words, there’s only one group layer per directory and the hierarchy of the directory isn’t preserve.
The expected result is shown in the images below.

Here's my code:

import arcpy
import os

folder =r'C:' #root folder
template_group_layer = r'C:' #empty layer
mxd=arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

groups={}
for path, dirs, files in arcpy.da.Walk(folder): #search in directory and create group layer

for d in dirs:
    lyr=arcpy.mapping.Layer(template_group_layer) #add empty group layer to ArcMap
    lyr.name=d
    arcpy.mapping.AddLayer(df,lyr)
    groups[d]=arcpy.mapping.ListLayers(mxd,d,df)[0]

for f in files: #select files that are shapefiles and insert them in their groups as the directory architecture
    fp=os.path.join(path,f)
    dsc=arcpy.Describe(fp)
    print dsc
    lyr=None
    view=None
    if dsc.dataType == 'ShapeFile':
        lyr=arcpy.management.MakeFeatureLayer(fp,os.path.basename(fp))[0] #edit 1 - remove layer in dataframe to prevent duplicate
        arcpy.mapping.RemoveLayer(df, lyr)
        print lyr
    else:continue

    if path==folder and lyr:
        lyr.visible=False
        arcpy.mapping.AddLayer(df,lyr)
    else:
        d=os.path.basename(path)
        arcpy.mapping.AddLayerToGroup(df, groups[d], lyr, "BOTTOM")

Current result


texts: Tenant isolation in SaaS product based on PostGIS - do the same concepts as PostgreSQL apply?

texts:
We're working on two different GIS apps and we're a bit torn around how to provide a level of tenant data isolation on our GIS systems.

System A would need to have strong tenant data isolation where 1 single app is used by multiple tenants to access their own separate data. We'd be looking at 50-500 range of clients to support and therefore up to 500 individual databases to connect to.

Our plan so far is to find a way to work with isolated shards or altogether different databases. Do all the capabilities of a standard PostgreSQL apply in this scenario? if anything works with PostgreSQL are we safe to assume it will work with PostGIS?

System B does not need complete tenant data isolation but it does need query-level isolation whilst keeping it possible for a 'tenant' manager to query all data.

This is aimed for a single company which has 50+ projects with separate users that need to query only project-relevant data, whilst at the top level being able to also query the data as a whole. We'd like to optimize for project-level query performance in this scenario.

So there's a project-level isolation when running queries for business/performance and to some degree security purposes.

Has anyone encountered any of these use cases with PostGIS/GIS data? We've yet to make a choice for database but PostGIS seems the obvious choice so far.


texts: How do I make gdalwarp set target extents to -180 -90 180 90?


texts: How to perform intersect operations in the Field Calculator in QGIS?

texts:
How to use the Intersect function?

How would I specify to

REPLACE layer1.field1 with layer2.field2 WHEN layer1.geometry INTERSECTS layer2.geometry

Is that feasible in QGIS?


texts: OpenLayers serving a MapServer layer: data won't display correctly without using WMS

texts:
I am trying to serve my Mapfile data as a MapServer layer in OL, rather than a WMS layer.

When I display it using WMS, everything works fine. The part of the mapfile I have struggled with was SIZEUNITS. Seems easy enough, but I encountered some strange behavior along the way.

The issue: my SIZEUNITS units is miles and the size is 100 (so a 50 mile radius). When I use the WMS layer below, everything is fine:

var cust = new OpenLayers.Layer.WMS("POIs",
"http://192.168.2.164/cgi-bin/mapserv.exe?map=/data2/html/html2/david_POIs4.map",
{layers: ['POIs']}, {wrapDateLine: true, singleTile: true});

However, when I try to switch it over to:

var cust = new OpenLayers.Layer.MapServer("POIs",
"http://192.168.2.164/cgi-bin/mapserv.exe",
{map: "/data2/html/html2/david_POIs4.map",
layers: ['POIs']}, {singleTile: true} );

The points (circles rather, with points in the middle) are HUGE. More or less taking up the whole map.

So, why would switching from WMS to MapServer show such glaring differences? A MapServer layer should be simpler to correctly configure.

When I was stuck on a similar issue last week I posted here (of course hadn't gotten this far yet).

Edit: see the images here for a visual. First map is correct (zoomed in), 2nd map is incorrect.

Edit2: I've never had an issue with projections before (I have the same projection set in both the LAYER and the MAP). Here's the layer in question, and if I use SIZEUNITS MILES it seems to automatically get blown out of proportion, regardless of the value you pass to SIZE within STYLE.

    LAYER
            NAME POIs
            DEBUG 5
            TYPE POINT
            CONNECTIONTYPE OGR
            CONNECTION 'PointsOfInterest_Latest.db'
            DATA 'POIs'
            STATUS ON
            SIZEUNITS miles
            LABELITEM "name"
            CLASS
                    STYLE
                            SYMBOL 'Airport'
                            COLOR 8 29 118
                            #WIDTH 5
                            SIZE 10
                            #MINWIDTH 1
                            #MAXWIDTH 1
                    END
                    #STYLE
                    #       SYMBOL 'City'
                    #       COLOR 8 29 118
                    #       SIZE 6
                    #END
                    LABEL
                            FONT "andalewt"
                            TYPE truetype
                            ENCODING "UTF-8"
                            MINSIZE 9
                            MAXSIZE 9
                            POSITION AUTO
                            PARTIALS FALSE
                            COLOR 8 29 118
                            OFFSET 1 1
                            PRIORITY 10
                            BUFFER 0
                    END # LABEL
            END

            PROJECTION
                    "proj=latlong"
                    "ellps=WGS84"
                    "datum=WGS84"
            END
    END

Now, you can see the SIZEUNITS set to MILES and the SIZE set to 10. Being the farthest zoom level out, you should barely be able to see the circles, but they are huge, as you can see here.

Any input is appreciated.


texts:
I found a sample of code

SELECT *
FROM locator
WHERE SQRT(POW(X(center) - 49.843317 , 2) + POW(Y(center) - 24.026642, 2)) * 100 < radius

in this question Find points within a distance using MySQL.

How can I add an ORDER BY distance to this code?


texts:
I am working with ArcMap 10.3.  I have imported DEM files.  When I imported the files they had the extension of .DDF.

I would now like to be able to take a list of XY coordinates and, for lack of a better term, merge them with the DEM files so that I can get the elevation of any given point.

I imported XY Data and projected it with the NAD 1983 coordinate system.  But, I can't figure out how to "combine" the 2 layers to get the elevation of my XY points.  I can right click on the map and click "Identify..." and get the elevation, but that's about it.

Also, if there is a way to do this without the Spatial Analyst extension, that would be great.  I am supposed to have it on my license, but for whatever reason I don't and need to know if there is a way to do this without it.


texts: Python 3.6, Django 1.11, Windows 10, Gdal 2.2.4 not working

texts:
I've been trying to find a solution to this issue, but couldn't find anything.
I've installed OS4GEOW using the Django Docs

But I get the following error:

  File "C:\Python36\lib\site-packages\django\contrib\gis\gdal\libgdal.py", line 49, in <module>
    lgdal = CDLL(lib_path)
  File "C:\Python36\lib\ctypes\__init__.py", line 348, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: [WinError 127] The specified procedure could not be found

I've modified the C:\Python36\lib\site-packages\django\contrib\gis\gdal\libgdal.py file, so it contains this line:
lib_names = [str('gdal201'), str('gdal202'), str('gdal111'), str('gdal110'), str('gdal19')]

instead of:

lib_names = [str('gdal201'), str('gdal20'), str('gdal111'), str('gdal110'), str('gdal19')]

gdalinfo --version returns: GDAL 2.2.4, released 2018/03/19

I also have GDAL=="2.2.4" python package installed.

My path contains:

C:\Python36\Scripts\
C:\Python36\
C:\OSGeo4W\bin

My C:\OSGeo4W\bin folder contains gdal202.dll which in details says 2.2.4 version.

But it still doesn't run, would be glad for anykind of pointers/help because I'm kind of lost with GDAL Python 3.6 combo.


texts: FDOGraphicsLayer IGraphicsContainer Reset/Next throw 'not implemented'

texts:
Running ArcGIS 9.3.2 w/ ArcFM.

protected override void InternalClick()
{
    var d = _App.Document as IMxDocument;
    var m = d.FocusMap;
    var c = m as IGraphicsContainer;

    c.Reset();
    for( var e = c.Next(); e != null; e = c.Next() ) {
        var p = e as IElementProperties3;
    }
}

If the Active Annotation Target is set to anything other than '<Default>' then c.Reset() and c.Next() throw 'not implemented' exceptions.  I interpret this as FDOGraphicsLayer not implementing those methods.  Is this expected ?  I could not find any documentation that mentions this, nor does it seem correct.

I currently work around the issue by calling c.LocateElementsByEnvelope and pass it the ActiveView extents (basically locate all elements).  Casting m to IGraphicsContainerSelect and calling SelectAllElements also works.  Both seem like hacks.


texts: Google Earth Engine : error: Cannot use 'in' operator to search for 'type' in {"type":"LinearRing","coordinates":

texts:
I'm trying to export Sentinel-1 bands with this code in Google Earth Engine but I have recently this error:

  Cannot use 'in' operator to search for
  'type' in {"type":"LinearRing","coordinates":

I didn't have this problem in the past. So here is the code.

var start = new Date("02/05/2015");
var end = new Date("03/05/2015");

  var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
      .filterDate(start, end)
      //.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
      //.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
      // Filter to get images collected in interferometric wide swath mode.
      //.filter(ee.Filter.eq('instrumentMode', 'IW'))
      .filterBounds(geometry);

var count = s1.size().getInfo();

for (var i = 0; i < count ; i++) {
  var img = ee.Image(s1.toList(1, i).get(0));
  var geom = img.geometry().getInfo();
  Export.image(img.toFloat(), img.get('system:index').getInfo(), {
    'scale': 20,
    'crs': 'EPSG:32632',
    'region': geometry.toGeoJSONString()
});
}

Export.table(s1,'Sentinel1_csv');

print(Chart.image.seriesByRegion(s1, geometry, ee.Reducer.mean(),'VH')
.setOptions({title: 'VH mean value over time on R3'}));