texts: Image classification using google satellite images

texts:
Currently our expert user classifies locations by viewing the Google earth satellite image of a location. They consider circular areas of a 5km radius at a set zoom level. They  evaluate the image for density of buildings or trees in within the circle and assign the location a classification. We get data in the form:

location         | Classification
(lat-long)       |
----------------------------------
123.123, 12.34   | B
122.123, 12.34   | A
121.123, 12.34   | C
120.123, 12.34   | A
119.123, 12.34   | A

I want to use this as training data to build a system that can automatically classify a location based on the satellite image.

To train the system, I would need to extract relevant parameters from the satellite images.  I guess the image parameters would need to be rotation independent. Possible parameters:

total intensity for RGB values,
histogram
count of edges
??

I plan on using C# and available AI & image processing libraries (e.g. opencv, encog, AForge) to do implement the training and classification stages.

Some of the image processing tutorials concentrate on object recognition but that's not really what I'm after (I think). It's more about grouping the images into different categories based on more general properties of the image.

Can anyone advise on a) the feasibility of this plan and b) what parameters I should be extracting from the images?

Thanks


texts:
I am new to arcpy, and I am trying to write a code in order to extract just the records containing specific attributes.
This code gives back an empty table. What am I doing wrong?

import arcpy
arcpy.env.workspace = 'C:/work'
arcpy.env.overwriteOutput = True
fc = 'rbh.shp'
exclude = ["02-R00-04-B$0$210", "02-R16-04-B$0$210", "02-R00-04-B$0$220", "02-R16-04-B$0$220", "02-R00-04-B$0$211", "02-R16-04-B$0$211", "02-R00-04-B$0$320", "02-R16-04-B$0$320"]
with arcpy.UpdateCursor(fc, "Layer") as rows:
    for row in rows:
        if row not in exclude:
            rows.deleteRow()


texts: Combining rasters to calculate total inputs

texts:
I am working on creating a map of total chemical inputs into a study area.  I have several raster layers of different inputs with values stored in the attribute tables.  I am having a hard time figuring out how to add these inputs together.  Is there a way in the raster calculator to add a specific field from the attribute table of each raster (these values are all decimals) and to ultimately display total chemical inputs spatially?

** I am using ArcGIS 10.3 with spatial analyst. **


texts: How to make sure labels in qgis turn up in layout map?

texts:
I am using natural earth for qgis 3.0. I have rule based layers for 10m_admin_0m_maps_unit to show only a certain number of countries in my Caribbean map. That now looks perfect when I am in the qgis, but when I switch to the layout view to build my map, the layers are gone.

I have looked for some tutorials and saw that I should do the following
click layer properties, then under rendering, select show all labels for this layer (including colliding labels).
I have done this and refreshed my layout (the map is not locked), but the labels still do not appear.

I see that I can manually add the country names, but is there a way to get the selected names of this layer to show up on the layout map?


texts: Reducing time for Export to PDF (or any other alternative to print maps) in ArcPy?

texts:
For my work, I designed a tool to create PDF maps with large amount of data, intersected with 5 buffer rings using ArcPy. The job involves intersecting points, selecting points, generating Elevation for the points and distance from the order location (the centre of the buffers). The job works fine for single map and takes only 2-3 minutes.

But some other requirements came up and now I am appending two more PDF to the original PDF but with different extents based on the ring radii. The job is taking as long as 8-9 minutes. If I create mapbook with data driven pages, it is not that slow. I have tried different resolution, it is somewhat efficient with less dpi but project manager is not happy about it. We are using 400 dpi instead of default 300 dpi. Also, Layers and attributes both are required for reports- though that change did not effect time that much. The scaling of tool is unknown too. The original design was meant for 1000 times uses in one month, but business is launching in another jurisdiction, so not sure about the number of uses down the road. And with more uses and longer run time- the tool will be running into 'Execution of job XXX failed because of crash or termination of the server object.' regardless of ArcGIS service time change.

Has anyone used another PDF generation tool with in Python to generate maps?

if (country.lower()== 'xy'):
            scale = df.scale * 1.1
            df.scale = ((int(scale)/100)+1)*100
            MainTitleTextElement.text= "Map : "+ Buffer4 + " Mile Radius"
            outputLayoutPDF4 = os.path.join(scratch, 'mapExport4.pdf')
            arcpy.mapping.ExportToPDF(mxd, outputLayoutPDF4, "PAGE_LAYOUT", 640, 480, 400, "FASTER", "RGB", True, "ADAPTIVE", "RASTERIZE_BITMAP", False, True, "LAYERS_AND_ATTRIBUTES", True, 90)
            scale = df.scale * 1.1
            df.scale = ((int(scale)/100)+1)*100
            MainTitleTextElement.text= "Map : "+ Buffer3 + " Mile Radius"
            outputLayoutPDF3 = os.path.join(scratch, 'mapExport3.pdf')
            arcpy.mapping.ExportToPDF(mxd, outputLayoutPDF3, "PAGE_LAYOUT", 640, 480, 400, "FASTER", "RGB", True, "ADAPTIVE", "RASTERIZE_BITMAP", False, True, "LAYERS_AND_ATTRIBUTES", True, 90)
            outputPDF.appendPages(outputLayoutPDF4)
            outputPDF.appendPages(outputLayoutPDF3)


texts:
I have been needing to import anywhere between a few to 10's of text layers into a project per day for my job.

The text files are in a .log format (Space Delimited)
The X, Y, and CRS designations are the same for all of the files

Is there a plug-in or python command that would help me import these files more efficiently?


texts: Making land-use raster reclassification in QGIS?

texts:
I have been struggling with the task of creating a least-cost path between two points on a land-use map.
So far I have got the raster as follows (category names don't matter right now)

So these 7 categories are the values of my raster cells, right? (0-6). I would like to change these for:
0 to 100 (highest cost)
1 to 0 (lowest cost)
2 to 80 (high cost)
3 to 30
And so on... (the logic behind this would be of an animal crossing it, that's why it looks odd).

The problem is that I can't seem to find a way to reclassify this. I see people talking about this with DEM files, but I haven't found much information about how can I do this with a "simple" land-use raster. Am I trying to work with too little information?
I believe once I get this reclassified, it will be easy to calculate the least-cost path between two points in the raster, right?


texts: How to generate shapefile from Excel or CSV files?

texts:
Possible Duplicate:
  How can I convert an excel file with x, y columns to a shapefile?

Probably this could be an easy question for most of the people. But I would like to know the way to generate a shapefile, if I have an input as Excel or CSV file format.


texts:
I have a single 2.5 GB JP2 image of a city that I need to split into tiles, which I'll then process to create flight simulation scenery.

Google as I might, I can't find the solution.

The JP2 apparently has GIS XML data embedded, which doesn't affect me either way, since I don't need it.  All I'm trying to do is break the image down into x-number of preferably 2048 x 2048 tiles.

A utility to do this in either OSX or Win7 would be great.

Ideas?

Thanks!


texts: Get Raster Values from a Polygon Overlay in ArcGIS

texts:
I have the same problem as Curlew in Get Raster Values from a Polygon Overlay in Opensource GIS Solutions and thecrashlandingdodo in Step-by-step: How do I extract Raster values from Polygon overlay with Q-GIS or R?. Curlew says it's rather easily done in ArcGIS, but I just don't know how and R is on my to-learn list.

I have many polygons in a shapefile and a raster with different values. I would like to know for each polygon how many percent each raster value has in that polygon. Any of the zonal stats does not seem to do this (or did I miss something?).


texts: Remap with number-limits and not individual values in GEE

texts:
I have a rather basic question regarding the reclassification of ImageCollections in GEE. When I want to reclass individual values, I use the remap function as follow

var my_reclass= function(image) {
  var myreclass = image.remap([10 ,12],[  0,   3], null, 'my reclass');
};

What I would like to do though and I can not find any examples, is to use number limits for my reclassification

e.g.

-1 thru 3 = 5

3.1 thru 8.2 = 1

...

How can I achieve this?


texts: Using Intersect Tool in ArcGIS gives empty output?

texts:
I'm currently searching for correlations between pavement conditions and crashes on highways.

I created a 420ft buffer layer for my selected highways and now I am intersecting it with  a layer for crash points. By looking at the map I can clearly see that there are points which intersect with my buffer layer, but The new intersect layer is empty with nothing in the attribute table.

I am getting a warning-000117 from Intersect.

Did I do something wrong?

I'm certain I followed the correct procedure for buffering and intersecting, but I end up with nothing.


texts: Importing an image from hyperlink to my table of contents (layer) in ArcMap

texts:
I have been trying to find a way to make it possible to import a file from my computer to the table of contents via a hyper link.
For example: We have organized a grid, so when I identify I get information of the file path to the Orthophoto for that grid square. When I manually have to go and find that folder and file among thousands of others. So instead of having to find the file and import it I would like to identify the grid square and that grid square should have a hyperlink as a field that, when clicked adds corresponding file to my table of contents.
I have quite a limited knowledge of scripting, but I have spent hours of googling and have not found any way to do it, the only thing I find is examples of scripts for opening files in web-browser and in external softwares as windows image viewer.


texts: Converting *.asc to *.shp without using ArcGIS for Desktop or QGIS?

texts:
I am trying to perform file conversion programmatically, using C#.

I have read the ESRI Shapefile Technical Description (An ESRI White Paper—July 1998).

I cannot calculate records for *.shp file from *.asc file.


texts: Find area under different colors in TIFF?

texts:
I have a TIFF file that has 17 bands (colors). Each color corresponds to different land use land cover characteristics.

Additionally, I have a shapefile of a district with blocks. I want to find the area covered under all 17 colors in all blocks. I can use zonal statistics or r.report. But both of them work for a single district shapefile not for district shapefile with blocks. What can I try?


texts: UpdateCursor if loop references next row in field

texts:
I am trying to populate a unique identifier field "STD_ID" using arcpy.da.UpdateCursor.  The identifier would be an iterative 1,2,3 and would reset the iteration based on region ID in which the polygon is located "RBSB_ID2"

My problem is I can't figure out how to reference the next row to restart the iteration. I updated the for loop to include cursor.next(), the error is now that I never started the iteration.

#call stands layer and define variables
stands = "F:/ArcGIS/Data/ForestPlanning.gdb/NYC_Stands"
rbsb_std = ('RBSB_ID2', 'STD_ID')

#define iteration.
iteration = 0

#create cursor for Stand IDs using arcpy.da.UpdateCursor

with arcpy.da.UpdateCursor(stands, rbsb_std) as cursor:
for row in cursor:
    try:
        next_row = cursor.next()
    except StopIteration:
        next_row = row
    if row[0] != next_row[0]:
        row[1] = 1
    else:
        row[1] = iteration + 1
        iteration = iteration + 1
    cursor.updateRow(row)

Any ideas on how to fix this?

I just edited the code to reflect toslunds original answer.

In response to the desired output ... assume

RBSB_ID2 = [100023, 100023, 100023, 100144, 100144, 1000144, 1000144, 1001512, 1001512]

#The output for 'STD_ID' would be:

STD_ID = [1, 2, 3, 1, 2, 3, 4, 1, 2]

#The next step would be to convert to STD_ID to a string & later convert to something like

#[001, 002, 003, 001 ...] and then later concatenate with the RBSB_ID to
#create a unique id for each row that links RBSB_ID2 with STD_ID.

#RBSB_STD_ID = [100023001, 100023002, 100023003, 100144001, ...]


texts: GDAL.grid creates an empty tif file (with no error during computation)

texts:
I'm attempting to create a raster file from XYZ data. I used GDAL.Grid on a .csv file with XYZ data and .vrt definition. The code produces a .tif file with no data. I'm having a hard time troubleshooting my code as Python GDAL is not throwing any error on Google Colab. I tried changing data, .vrt definition, and parameters of gdal.grid.
A reproducible example is located here:


texts:
I have a very large GeoJSON file and I would like to make it faster.
How can I add in the filter a check to see if points are inside the map's view box, so only points in the current view are rendered?
 var str_pol = L.geoJson(null);
// load GeoJSON from an external file
$.getJSON('https://url', function(data) {
  var strutt_poly = L.geoJson(data, {
      style: subteStyle,
      onEachFeature: function (feature, layer) {
        var popupContent = "
Desc: "+ feature.properties.desc+'
';
            layer.bindPopup(popupContent);
      layer.addTo(strutture_pol);
    }
  });
});

I tried to use the WFS service.
I have GeoServer to another server I can see the WMS layer but not WFS layer, why?
    var sost_wms= L.tileLayer.wms("http://url:8080/geoserver/test/wms", {
        layers: 'test:sostegno',
        format: 'image/png',
        transparent: true
    });

var case = L.geoJson(null);
// load GeoJSON from an external file
$.getJSON('http://url:8080/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=test%3Acase&maxFeatures=50&outputFormat=application%2Fjson', function(data) {
  var casa= L.geoJson(data, {
      style: function(feature) {
      return {
       weight: 1,
                opacity: 1,
                color: '#000',
                dashArray: '3',
                fillOpacity: 0.5,
                fillColor: '#ec7f1a'
      };
    },
    onEachFeature: function (feature, layer) {
        var popupContent = "
nome: "+ feature.properties.nome+'
';

            layer.bindPopup(popupContent);
// aggiungi i dati del geojson al layer var clusters
      layer.addTo(case);
    }
  });
});

On my local server, I can see them both.


texts:
How do we make it possible to click, drag, copy text that appears in Info panel in leaflet.js?

In a leaflet application similar to the choropleth example, I'm populating the info panel on click instead of on mouseover. I need to enable the user to select and copy text that is present on the panel. At present, hyperlinks I make on the panel are clickable. But the text in the panel cannot be selected.. the mouse behaviour is like it's on the map only.. the panel is a "ghost" for the mouse unless a hyperlink comes along.
Link to project.


texts: Why can't Carto recognize the coordinates of my csv file?

texts: Carto can't recognize the point geometry of some of my csv files. For example I can visualize in ArcGIS the point geometry of this csv file I have posted here. But when I upload the file to Carto it doesn't work it can only recognize the attribute data in the table but not the coordinates in the XY coloumns. I have another file where it works and the XY columns look identically which puzzles me. I have tried text and double as field type. And english and german comma as well ( "," and ".") All four versions don't work. What could be the reason?


texts: What tool would you use to "clip" a vector line file to a vector polygon file?

texts: I am trying to select the use of a large vector line file within only the parameter of a vector polygon file to do spatial analysis. It would be the same thing as using the extract> clip tool under the analysis tools that is for two polygon files. What tool do I use that works to clip a vector line layer to a polygon layer?


texts:
I created a QGIS 3 Plugin with the Plugin Builder, on Linux, and wrote some code for it, and it works fine so far.

Now, to improve the readability and structure, I want to split the code up, in some modules. But if I then want to start QGIS and check if it works, the module isn't found. It's in the same folder than the other files and I tried to import it like

import module
from module import f1
from .module import f1

Where f1 is the only function I have defined so far, and contains only a print.
Also I'm working with Make and use  make deploy

EDIT:
If I import it in the IDE it works and I can use the function f1.
Also if I build it as a package with an __init__.py and then import it it works... just if I want to test it in QGIS it doesn't find the files.

SOLUTION:
Check your Makefile... configured my Makefile wrong, so the modules wouldn't get deployed to the QGIS plugin folder, hence it worked in my IDE but not in QGIS...
Thank you for your Ideas and Help!


texts:
I have model plotted in R

base<- raster("example.tif")
 plot(base)

image
I'd like to change value (0,5,10,15,20,25)

breakpoints<-c("0,5,10,15,20,25")
 plot(base,break=breakpoints)

Error in .rasterImagePlot(object, col = col, xlab = xlab, ylab = ylab,  :
  object 'breakpoints' not found

I'd like make ramp(default colour)


texts: Plotting PostgreSQL PostGIS points on map for web application

texts:
My goal is to create a web app where there's a map of a city and you can filter for year, month, day, etc. and crimes for that day will pop up as points on a map. I already have the data formatted as PostGIS objects in PostgreSQL, with the associated crime type and date, but I'm not sure how to go about creating the web app and querying the data on the front end.
I know that Python and possibly GeoDjango could be an option. I've also looked into gmplot. I want something that's dynamic and users could zoom in and click on individual data points to get some more information about the crime.


texts: Overlay two rasters and extract values from one of the rasters

texts:
I'm creating a Higuchi viewshed and I'm not sure what tool to use to do the last step.

Stages in the Generation of an Higuchi Viewshed

Calculate binary viewshed for viewpoint
Calculate a distance layer from the viewpoint
Reclassify the distance layer according to the Higuchi class criteria (determined by the typical tree height for the area and period under study)
Overlay the binary viewshed upon the distance layer to extract the distance range values which fall within the in-view zone

I've created the binary viewshed and the distance layer (using euclidean distance tool), and reclassified the distance layer.

However, I'm not sure what tool to use to extract the values from the distance layer that fall within the viewshed zone.

Example of the process:


texts: Creating multiple non-overlapping rings but unable to tell which fc they surround using ArcPy?

texts:
I am using ArcGIS Desktop 10.4.

I would like to create multiple ring buffers surrounding my polygon localities. My buffers should not overlap and I need to know, which location each of my buffers surrounds.

I've thought to simply use arcpy.MultipleRingBuffer_analysis. However, if I use Dissolve_Option = "NONE", I am able to define which locations my buffers surround, as the locations ID is stored in newly created "ORIG_FID" field. However, my buffers overlap.

If I use Dissolve_Option = "ALL", my buffers do not overlap. But, I have no information about which location does my buffer surround.

How do I accomplish having multiple distances non-overlapping buffers but be able to know the location they belong to?

My code to create multiple ring buffers:

# Import modules
import arcpy, os

# Set environmental settings
inWD = "C:/Users/Localities"
arcpy.env.workspace = os.path.join(inWD, "analyzed.gdb")
outPath = os.path.join(inWD, "output.gdb")

# Allow files to overwrite
arcpy.env.overwriteOutput = True

# Define input variables
inFC = "NPR"

try:

    # Define Multiple ring buffer variables
    distances = [100, 500, 1000, 2000]
    buffUnit = "Meters"
    Field_Name = "distance"
    Dissolve_Option = "ALL"
    Outside_Polygons_Only = "OUTSIDE_ONLY"

    # Process: Multiple ring buffer
    outBuff = os.path.join(outPath, "buff_complete")
    arcpy.MultipleRingBuffer_analysis(inFC, outBuff, distances, buffUnit,
                                      Field_Name, Dissolve_Option,
                                      Outside_Polygons_Only)
except:

    print(arcpy.GetMessages(0))

Simple illustration for distances [100,500,1000,2000] Meters: even if I click only on 100 m buffer (the closest to inFC), my buffers have to overlap, as in attribute table I've selected 3 different buffers.


texts: Where can I get a WMS or WFS for near-real time Wind (atmospheric) data?

texts: I've been looking for a good operational WMS or WFS free service to access near-real time atmospheric data in QGIS, namely wind vectors. I haven't found anything good yet.


texts: How to select features based on multiple values in graphical modeler?

texts: I'm trying to create a basic model using QGIS, the QGIS graphical modeler and the sextante toolbar. In part of the model, I need to select features based upon multiple attribute values (for example - A = 'value' OR A = 'another value'). The select by attribute algorithm provided by QGIS only allows you to enter a single value at a time - is there a simple tool/method to enter an expression to select multiple values from a dataset, and which can be incorporated in the graphical modeler?


texts: Interpolating irregular shape, merging DEM's using ArcGIS Desktop

texts:
My problem is quite complex, but I hope that solution may be simple or at least possible.

I do have .img DEM 0.5x0.5m. DEM consists of river bed, river valley and adjacent areas (it is strip of 5km with the river in the centre).
DEM was derived from ALS, so it does not show depth of water feartures and of course that something I am interested in. I do have some cross-section of the river bed and will have surveys of certain sections, lets say I have various point data for cross-section of the river bed, which will be 200-300m apart.
My aim is to interpolate points for those 200-300 m gaps, convert it in to DEM and somehow combine with the original DEM. So as you see there are several issues.

How to extract shape of the river corridor from original DEM?

How to interpolate points within such irregular shape as river corridor, so I can receive valid DEM for river bed?

How to combine original and river bed DEM?


texts: Applying rules to FrameCloud popups in OpenLayers

texts:
I have style rules for features that enable the map to display the features based on its rule. For example, the map displays a feature from specific scale.
I do this in this way

newRule = new OpenLayers.Rule({
    name: 'rulename',
    minScaleDenominator: 2000,
    maxScaleDenominator: 19000,
    filter: someFilter,
    symbolizer: someSymbolizer
});

Then add the rule to style and apply the style to features. This lets the map to display the feature from scale 2000 to 19000 and when the map goes out of the set scale then it does not display the features.

What I want to do here is, to apply those rules to FrameCloud Popups.

From what I've found so far, popups cannot have style rules. So I implemented this by registering event 'zoomend' and let it check whenever zoom level changes so that at certain zoom level the map loops through all the popups and hide / show them.

It works just enough but it seems unreliable and hardcoding.

I think there should be a better way to achieve this. Please advice me if you know the answer. Thanks in advance :)


texts: How to transform this PostGIS type of projection to EPSG:900913 in OpenLayers?

texts:
I have transferred this kind of projection from this type(which I don't what type it is):

010100000054E077E6F8F33640FD87F4DBD7AD4340

to EPSG:900913. I have the points which are in the above format in OpenLayers I just don't know the projection type so I can fidn the right command to change it.

So I am asking for both the type of the projection and any commands that will change it to ESPG:9000913 and add it to my layer.


texts: I've searched all over google and have found nothing. I want some tiff files I have in ArcMap to be in greyscale with the rest of the layers in colour, similar to the style override > greyscale option in MapInfo. I will be printing A1 posters and do not want the background imagery in colour. I've tried changing the symbology to stretched (as suggested on a forum), but it looks weird, like this:


texts: QGIS: Prevent the Primary Key from incrementing when feature won’t be saved

texts:
So, I have a PostgreSQL/PostGIS DB attached to my QGIS project (MacOS, 3.10.).

I start adding a feature and cancel the process w/o saving the form. When adding another feature my primary key increments by one even though the last feature wasn't saved.

Is this avoidable? Meaning: Can I prevent the PK from being incremented when canceling a feature addition?


texts:
I'm developing a Combo Box for an ArcMap toolbar, that at present simply lists all of the layers that are in the table of contents. However the height of the drop-down menu is very small, and I cannot see a property that I can set to change the default height of the drop-down list. I've attached a screenshot to demonstrate the problem. You can see that only one item is visible in the drop-down menu, although there are four layers in the TOC. The other layers can be accessed (look at the tiny arrow buttons, which do work), but it is really difficult to do.

In the add-in config.xml there is a rows property, which defaults to 4. Changing it to a higher or smaller number does nothing.

In the Combo Box documentation, there are properties for the Python class controlling the Combo Box. These allow control of dropdownWidth and width, bizarrely controlled by supplying a string of lots of the same character.

I've been using the Python add-in wizard to form the template. I'd like to stick to that but I'm open to ArcObjects if that's the only solution.


texts: Combining multiple netcdf in folder in R

texts:
I have a folder that contains 61 netcdf files containing 4 dimension (lat, lon, time and depth) and 9 variables.

I would like to combine all the files in the folder to one netcdf file with the variable "thetao" so that the 4 other dimensions remain.

I do not have a lot of background in R. So far I have come to this from what I have gathered from other sources:

library(ncdf4)
library(raster)

# Set working directory
setwd("C:/Users/L1586/Downloads/CMEMS_2015/05/BAL-HBM_PHY-201505-DailyMeans/")

# Declare data frame
df = NULL

#Open all files
files = list.files("C:/Users/L1586/Downloads/CMEMS_2015/05/BAL-HBM_PHY-201505-DailyMeans/",pattern='*.nc',full.names=TRUE)

# Loop over files
for(i in seq_along(files)) {
  nc = nc_open(files[i])
}


texts: Run field calculator script in Python window

texts:
I have a script that runs perfectly in the "Field Calculator." The codeblock looks like this:

def timeday(time, am, diur, pm, noc):
    if(time >=am and time < diur):
        return "crep_am"
    elif(time>= diur and time<pm):
        return "diurnal"
    elif(time>=pm and time<noc):
        return "crep_pm"
    else:
        return "nocturnal"

and the expression looks like this:

timeday(!LMT_TIME!, !B_C_AM1!, !B_DIUR1!, !B_C_PM1!, !BEG_NOC1!)

As I mentioned, I ran this on a single attribute table for a shapefile and it worked perfectly, but I've never used the python window before. I want to run this code on all of 100 shapefiles. I've been playing around with it, but I can't seem to get it right. Here's what I have so far:

import arcpy
arcpy.env.workspace = r"C:\Users\evan\Desktop\Evan_4-
2017\Deer_Data\LaSal_Deer\GPS_Data\GEO_WGS1984"
Shapefiles = arcpy.ListFeatureClasses("junktest.shp", "junktest2.shp")
codeblock= """def timeday(time, am, diur, pm, noc):
if(time >=am and time < diur):
    return "crep_am"
elif(time>= diur and time<pm):
    return "diurnal"
elif(time>=pm and time<noc):
    return "crep_pm"
else:
    return "nocturnal""""

for shp in Shapefiles:
    arcpy.CalculateField_management(shp,"TIME_DAY","timeday(!LMT_TIME!, !B_C_AM1!, !B_DIUR1!, !B_C_PM1!, !BEG_NOC1!)","PYTHON",  codeblock)

Alternatively, if there is a way to run the field calculator in ArcMap on all of the shapefiles as a batch, I'd be fine with that, but when I tried doing that I was not able to enter my codeblock, only the expression.


texts:
I have a series of layers of lines (shapefiles) in My Project. The CRS of some are different and to merge them they all must have same CRS.

When I use Processing/Toolbox/Qgis_algorithims/Reproject_layer the reprojection only works if it is allowed to save to temporary file somewhere (it will not change the CRS when I reproject if I try to save and replace it in My Project directory); THEN

I must remove the original from the layer list and reproject the temp file (without changing the CRS) to get it into my project and name it appropriately.

Is this the only, or proper way to go about this?


texts:
New to python scripting.  One of the issues I have been working on is trying to get a script to calculate the elevation gained in each stream section based on a single input and the knowledge of the elevation at that height. I have the code written so it can make all the correct selections but my issue comes from trying to get the script to fill in the attribute table after.  If anyone has used the book Python scripting for arcgis this is the exercise 2-8 in the book.
 Fields -
Elevation gained (this is a field that contains a single number about how much elevation is gained in each section.)
Flowline_out (this is the bottom section of every segment of creek.  The first one is a known number, after that the Flowline_out is the result of the last segments flowline_in)
Flowline_in ( this is the top part of each segment.  This number is generated from the addition of elevation gained to the flowline_out)

I am trying to use the update cursor

fields = ("Flowline_In","FLowline_Out","Elevation_Gained")
        elcursor = arcpy.da.UpdateCursor("AnalysisCreeks",fields)
            for row in elcursor:
                i = arcpy.getValue("Flowline_Out") + arcpy.getValue("Elevation_Gained")
                row[0] = i
                cursor.UpdateRow(row)
            del elcursor, row


texts: Automatic selection of legend object in QGIS

texts:
I'am using the Corine Land Cover data for a projet. There is around 40 types of land use in this database. In my print composer with my scale (1:20 000) I neve have more than 10 type of land use.

I would like to show in the legend, only those 10 type of land use. Is that possible? (This option is available for the table)

I'am using QGIS master 1.9.


texts: Detecting self-intersecting polygons in OpenLayers

texts:
I am using OpenLayers Version: 6.3.1. I am allowing users to draw a Polygon in my website. I want to validate the polygon by checking for self intersection. Like This:

I referred to Detecting polygon self-intersection and proceeded with using turf.

 draw.on('drawend', function(){
      var layerSource = vector.getSource();
      layerSource.forEachFeature(function(feature){
        var poly = feature.getGeometry().getCoordinates();
        var kinkedPoly=turf.polygon(poly);
        console.log('kinkedpoly='+kinkedPoly);
        var unkinkedPoly=turf.kinks(kinkedPoly);
        console.log('UnkinkedPoly=' + unkinkedPoly);
      });
});

I don't know how to proceed further. What should I change? Can someone suggest any other method?


texts: Windows cursor being outside of visible bounds of map using ArcObjects

texts:
What is the method to determine whether the Windows cursor is outside of the visible bounds (envelope) of the map?
I need to know if the cursor is currently over another window frame and not the map itself (i.e. the cursor is over the TOC, menus, or outside of ArcMap all together). I have used ArcMap.Document.ActiveView.ScreenDisplay.DisplayTransformation.VisibleBounds to get an IEnvelope and determine the visible bounds of the map, but now I need to know whether the cursor is outside of that or not.  I realize this must be a simple concept, but I have not been able to find a working answer.


texts: Are State Plane Coordinates always derived from longitude and latitude?

texts:
Do state plane coordinates always begin as longitude and latitude?

In other words, they are not another system that can be calculated without Long. and Lat., correct?

I guess I'm saying, back in 1930+ they came up with the State Plane system.  To make that work, however, they had to start with L&L.  Today, I am guessing I can simply compare it to a database.  But if for example, I have a GPS device and I want to know the Plane Coordinate, I have to convert it.


texts: GEE ：Cloud/Cloud Shadow Masking for Landsat TOA Data

texts:
Due to data limitations, I need to use TOA data to extract cloud shadows instead of SR data. Is there any good way?
In addition, why is the amount of SR data less than that of TOA?
and this is the code I used.
var qa = image.select('BQA');
var pattern = ee.Number(2).pow(4).toInt();
var mask = qa.bitwise_and(pattern).rightShift(4);
return image.updateMask(mask.not());


texts: Buffer view that calculates numbers from a different layer

texts:
Basically people draw lines and from these lines, I need to have a view that will create a buffer that will calculate potential revenue from those lines

mylinetable is a table with a portfolio_id and a line geom

mypointtable is a table with portfolio_id, units and a point geom

I've read and tried to adapt it for my code with this question but no luck so far How to create dissolved buffers in QGIS with PostgreSQL/PostGIS?

My current code (not working)

 CREATE or replace view myview AS
SELECT
a.path[1] as gid,
sum(b.units) filter (where b.portfolio_id = a.portfolio_id) as same_portfolio,
sum(b.units) as total,
a.portfolio_id,
a.geom::geometry(Polygon, 27700) as geom
FROM
(SELECT portfolio_id,
 (ST_Dump(ST_UNION(ST_Buffer(geom, 250)))).* as geom
FROM mylinetable group by portfolio_id
) as a left join mypointtable b on st_contains(a.geom, b.geom)
group by a.portfolio_id

The intended outcome is to sum all the units from each point from the line, and give people a visual output of the lines they are drawing on QGIS.

Both tables have a portfolio_id, I need a sum of units when portfolio_ids match with the lines and a different column of simple totals

This view needs to be loadable on QGIS that's why I'm adding the path for an id.the row number could also work

Am I missing something or..?


texts: Intentionally offset GPS coordinates using heading

texts:
Is there a known way of intentionally offsetting a latitude/longitude position by using a combination of the latitude/longitude and heading? I'm struggling to find any kind of answer or solution for this...

[Edit Partial Solution]

I came across this link which is what my code below is based on below which I think was a link within the Richard's link below.

As can be seen in the image below, I wanted to purposely introduce an offset which can be seen in the points highlighted in 'yellow' are offset by about 10 meters.

This does not hold true for some coordinates (i.e. I don't observe a consistent offset), which I am assuming depends on using the heading value, which I believe comes back to my original question, I believe use of the heading in a calculation is required.

I will continue my search and update if i I come across a solution for the next person.

#include <math.h>
#define EARTH_RADIUS (6378137.0)
void function(double lat_in, lon_in, double offset_mtr, double *lat_out, double *lon_out)
{
    double dLat = 0.0;
    double dLon = 0.0;
    if( (lat_out != NULL) && (lon_out != NULL) )
    {
        dLat = (offset_mtr/EARTH_RADIUS);
        dLon = (offset_mtr/(EARTH_RADIUS * cos(M_PI*lat_in/180)));

        (*lat_out) = lat_in + dLat * 180/M_PI;
        (*lon_out) = lon_in + dLon * 180/M_PI;
    }
}


texts: Does anybody know what coordinate system this could be?


texts: How to calculate a position based on nearby sensors with weights

texts:
I'm working on an indoor-positioning project where the users are carrying phones which can detect beacons located throughout the facility.

The phone app only knows the IDs of the beacons, and can detect the beacons' relative signal strength. I know the XY coordinates of each of the beacons. From this information I need to determine the user's position.

In the example below, the user's phone detects these beacons, with these relative strengths:

BeaconID   Signal Stregnth
    1            10
    2            3
    3            7

I know the user is located somewhere near those beacons - is it possible to determine their position based on this scant information?

Note - as far as I can tell the signal strength is relative and doesn't correlate to a specific distance. I've drawn the user midway between the 3 beacons but it's possible that they're located outside the triangle


texts: Getting ArcGIS Templates in Visual Studio 2010?

texts:
I have a work computer which has ArcGIS desktop installed first and then VS2010 (just the way IT does it) and then templates are not available, anyone know what I need to do to see the templates?

I am running this under Windows XP.


texts: Using ogr2ogr to reproject and rename into shapefile?

texts:
I'm trying to reproject, move from tab to shp and rename all the .tab contained in a folder tree using ogr2ogr.

I'm working with Windows 7, QGIS 1.8 installed with advanced OsGeo4w and gdal 1.9.2. I have no experience in batch.

Here's the code :

for /R %%f IN ('D:\3.BATCH\FAUNE\*.tab') do C:\Users\Samy\Desktop\SIGMA\15.webmapping\ms4w\tools\gdal-ogr\ogr2ogr.exe -s_srs EPSG:3727 -t_srs EPSG:2975 %%f "ESRI Shapefile" %%f_2975.shp %%f
pause

By calling this .bat into the cmd line, nothing append. When I simplify the code, I have in the best case "Unable to open" as Nik posted in the 3rd ressource below.

I wrote this code with the following questions:

Using ogr2ogr to convert all shapefiles in directory?
Is there any way to reproject multiple layers in QGIS?
How to reproject shapefile using ogr2ogr

According to resource 1, I translate the -f "MapInfo File" "%~dpnf.tab" "%f" into %%f "ESRI Shapefile" %%f_2975.shp %%f because as it is shape, I can't use dpnf.

According to resource 2, I wrote -s_srs EPSG:3727 -t_srs EPSG:2975 to reproject.

Finally I checked that the "Hide extensions for known file types" was off, according to ressource 3.

I don't know what to do now. There seems to be so many different ways to write this one line of code.


texts: Recursively find all GeoTIFF in current and subfolders and clip using GDAL

texts:
I have thousand of GeoTIFF file in subfolders, and I would like to crop all of them using shapefile and save in a Crop folder.
My working directory is:
/Data/2000/image_20000101.tif
/Data/2000/image_20000102.tif
/Data/2000/..
/Data/2001/image_20010101.tif
/Data/2001/image_20010102.tif
/Data/2001/..
/Data/2002/..
/Data/polygon.shp
/Data/Crop/

My script:
% for i in `find . -type f -name "*.tif"`; do gdalwarp --config GDALWARP_IGNORE_BAD_CUTLINE YES -srcnodata -999 -dstnodata NoData -cutline polygon.shp -crop_to_cutline $i ./Crop/crop"_"$i; done

Seems the looping and finding all geotiff in all subfolders is worked, but it failed to save the output.
ERROR 4: Attempt to create new tiff file `./Crop/crop_./2000/image_20000101.tif' failed: No such file or directory`

I am expecting the output inside Crop folder
/Data/Crop/crop_image_20000101.tif
/Data/Crop/crop_image_20000102.tif
...

How to solve above problem?


texts: How can I tell the accuracy of coordinates conversion by cs2cs?

texts:
For example, when I try this:

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=tmerc +lat_
0=22.31213333333334 +lon_0=114.1785555555556 +k=1 +x_0=836694.05 +y_0=819069.8 +
ellps=intl +towgs84=-162.619,-276.959,-161.764,0.067753,-2.24365,-1.15883,-1.094
25 +units=m +no_defs

When I put

  114.1689 22.4518

It outputs

  835447.18       834705.40 3.11

Can I tell from the result, the accuracy and possible numeric error of the result? Some meters or some percent?


texts: Problem with view of OpenLayers plugin map and Swiss national coordinate system in print composer

texts:
My projects are in the Swiss coordinate system (CH0193+ / LV95) and I use OpenLayers Plugin for the base map (Google streets).  In the normal window, this seems to work fine:

As soon as I open print composer and export the map as PDF, this happens:

As you can see, vectors and map are not congruent any more.

Do you have any idea how I can solve this?

I'd prefer a solution without changing the whole grid, because all the data I am given is this Swiss format.


texts:
I use QGIS 3.10 as GIS client and i would like to be able to load tiles served from GeoServer.
I've well installed the Vector Tiles extension in GeoServer, my layers are well served when I try to access http://localhost:8080/geoserver/gwc/service/tms/1.0.0/.
I've installed the extension Vector Tiles Reader in QGIS.
However, when I add my url in Vector Tiles (add layer) i've got an error The field 'vector_layers' is required but not found. This is invalid TileJSON


texts: Set number of columns and rows with PyQGIS

texts:
I know how to resample a raster using pyQGIS 3.0 and the gdal warptranslate function by imposing the pixel dimension of the output file, see below.
rasterLayer = QgsRasterLayer('/path/inputfile.tif', 'raster')

parameters = {
    'INPUT': rasterLayer,
    'SOURCE_SRS': "EPSG:32632",
    'TARGET_SRS': "EPSG:32632",
    'RESAMPLING': 1,
    'NODATA': 0,
    'TARGET_RESOLUTION': 30,
    'OPTIONS': "",
    'DATA_TYPE': 6,
    'TARGET_EXTENT': rasterLayer,
    'TARGET_EXTENT_CRS': "EPSG:32632",
    'MULTITREADING': True,
    'EXTRA': "",
    'OUTPUT': "/path/fileout.tif"}

processing.run('gdal:warpreproject', parameters)

Is there a way to resample a raster by imposing instead the number of columns and rows (pixel dimension can varies)?


texts: How can I draw georeferenced WGS84 grid on Android?

texts:
I'm planning to display a grid which is georeferenced and show the lat/long of the particular location depending on zoom level.

I want to display something similar to what MOBAC has as an option called 'enable wgs grid'. Could anyone help me get started? Or any idea how can I draw georeferenced grid?


texts: Clipping Feature Classes by Data View Extents with ArcPy?

texts:
Have a fully functional script about 90% of the time. It is that 10% that is really bothering me and I cannot seem to deduce a solution.

Script will take the user defined dataframe view extents and clip all layers within extent to a geodatabase. That 10% failure rate occurs when clipping an extent that does not contain one of the features (gate_line, fence_line, or wall_line).

My thinking would be to limit input for clipping to only layers that are found in the view extents, this is where I logically fail.

The current script:

import arcpy, os
from arcpy import env

#overwrite permissions
arcpy.env.overwriteOutput = True

#set map document to current
mxd = arcpy.mapping.MapDocument("CURRENT")
path = mxd.filePath
path = path.replace(os.path.basename(path),"")
arcpy.env.workspace = path

#creates new file geoGDB
arcpy.CreateFileGDB_management(path, "CurrentExtents.gdb")
gdb = path + "CurrentExtents.gdb"
arcpy.env.workspace = gdb

#set dataframe extents as a polygon
df = arcpy.mapping.ListDataFrames(mxd)[0]
dfAsFeature = arcpy.Polygon(arcpy.Array([df.extent.lowerLeft,
                                         df.extent.lowerRight,
                                         df.extent.upperRight,
                                         df.extent.upperLeft]), df.spatialReference)

lyrs = arcpy.mapping.ListLayers(mxd,"",df)

exportLayers = []
for lyr in lyrs:
    if lyr.isFeatureLayer == True:
        layerRplc = lyr.longName
        replace = [".", "\\"]
        result = arcpy.GetCount_management(lyr)
        count = int(result.getOutput(0))
#Only include non-blank features
        if count != 0:
            for item in replace:
                layerRplc = layerRplc.replace(item, "_")
            outFC = os.path.join(gdb,layerRplc)
            exportLayers.append(outFC)
            if arcpy.Exists(outFC):
                arcpy.Delete_management(outFC)
            arcpy.AddMessage("Getting our clip on...." + outFC)
            clip = arcpy.Clip_analysis(lyr, dfAsFeature, outFC)
#removes empty fc's from gdb
            clipRes = arcpy.GetCount_management(clip)
            clipCount = int(clipRes.getOutput(0))
            if clipCount == 0:
                arcpy.Delete_management(clip)
            del replace
    del lyr

del mxd
del df

Script is run in desktop from catalog with a saved mxd directory as the defined workspace.

Just started coding recently!


texts:
I've got a problem using QGIS. Always when I try to launch the programme with the desktop icon, nothing happens. I also checked the Task Manager and it doesn't show any process running regarding qGIS.

The target field shows this: "C:\Program Files\Quantum GIS Lisboa\bin\nircmd.exe" exec hide "C:\Program Files\Quantum GIS Lisboa\bin\qgis.bat"

I had Microsoft C++ 2008 already installed, so I launched the installer and selected to repair it. However, QGIS still isn't launching.

I'm using Vista 32-bit.

I've tried to launch it manually, but then I get an error message telling me, that qgis.bat couldn't be found. I don't have any .qgs files on my PC, so I can't test it.

I can't find qgis.bat in the bin folder. But there's a file there called 'qgis.bat.tmpl'.


texts: Documentation on wildcard syntax, such as for use in arcpy.ListFeatureClasses()?

texts:
I was wondering if anyone knew where I could find some documentation on proper syntax for use of wildcards in arcpy, such as with the ListFeatureClasses function?

I have been unable to locate any thus far.


texts:
I can connect to Oracle database adding and creating spatial layer.

In the same database there are non-spatial tables but I cannot view any of them when connected to the database.

I need to add these tables to QGIS to perform joins with spatial layer.

How do I add non-spatial Oracle table to QGIS?


texts: Modify raster values within QGIS raster style dialog

texts:
I'd like to use the QGIS Raster Style dialog to display modified raster values. Specifically, I'd like to see a gradient using the square root of the raster's value.

I understand it's possible to use the raster calculator to transform raster values, but - if possible - I'd like to skip that step and apply the mathematical formula during the styling, and avoid creating a modified copy of the original raster.

This type of operation is possible with vector files, using data defined values for colour/width/etc. Is it possible to do this with rasters?


texts: Finding area/multipolygon bounded by three other multipolygons

texts:
I am trying to get my head around PostGIS and QGIS.  I'm trying to find the area (in white) bounded by three multipolygons (in purple) - refer to picture.  Each multipolygon only touches the other two multipolygons; there are no intersections.

More info:  The purple areas are defined as multipolygons in the PostGIS table.  I'm interested in just the coordinates for the white area.  The coordinate reference system is WGS84.

Basically, the purple areas are maritime boundaries/EEZs which have been defined in the postgis table.  The white area refers to international waters/high seas pockets between those maritime boundaries.  I'm trying to identify the coordinates of the white areas so that I can create records for them in the same postgis table that has the EEZs.  I'm building a web-app for monitoring vessel movements near real-time between EEZs and high seas pockets.


texts: Methods for measuring line/with buffer overlap

texts:
I have about 250 routes and I want to know how many meters of the routes overlap with each other relative to the total meters traveled.

I have a buffer around each route so it is really how much of the buffers overlap if that helps.

I can provide sample data if that helps.


texts: Overlay Google Maps to OpenLayers in GeoServer

texts:
I am new to GeoServer and I want to overlay Google Maps to my OpenLayers in GeoServer, I don't know how to start, do suggest me how to begin.

I want know which file to access, file extension, and what language is being used to incorporate it in GeoServer?


texts: Differences beetween proj4 string and WKT

texts: What are the differences between proj4 strings and well-known text in terms of the ability to represent a projection, other than for WKT being more human readable?


texts:
I am using openlayers 3 to try and get vector data from my Geoserver 2.12.0. I am using ol.layer.vectorTile and ol.source.VectorTile to build my vector layer. The main goal is to pull polygon vector data from geoserver so the client can edit/save the polygon. I am using this as a reference to pull the Vector Tile data. I have completed everything on the tutorial except the "Create OpenLayers application" section. This is my openlayer's javascript:

var longitude = -97.587453;
var latitude = 35.471398;
var zoom = 11;
var format = 'application/x-protobuf;type=mapbox-vector';
var layer = 'test:poly';
var projection_epsg_no = '4326';

var layers = [
    new ol.layer.Tile({
        source: new ol.source.OSM()
    }),
new ol.layer.VectorTile({
        style:defaultFeature,
        source: new ol.source.VectorTile({
            tilePixelRatio: 1,
            tileGrid: ol.tilegrid.createXYZ({maxZoom: 19}),
            format: new ol.format.MVT(),
            url: 'http://192.168.10.15:8080/geoserver/gwc/service/tms/1.0.0/' + layer + '@EPSG%3A'+projection_epsg_no+'@pbf/{z}/{x}/{-y}.pbf'
        })
    })
];
var projection = "EPSG:4326";
var latLong = ol.proj.fromLonLat([longitude, latitude]);

var myControl = new ol.control.Control({target: $('#map-panel')});

var controls = [ol.control.Zoom, myControl];

var map = new ol.Map({
    layers: layers,
    target: 'map',
    control: controls,
    view: new ol.View({
        center: latLong,
        zoom: zoom
    })
});

This is the error I get from Chrome's Developer Tool Console:

  Failed to load http://192.168.10.15:8080/geoserver/gwc/service/tms/1.0.0/test:poly@EPSG%3A4326@pbf/11/472/1228.pbf: No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'http://192.168.10.15' is therefore not allowed access.

If I change the url to:

url: '/geoserver/gwc/service/tms/1.0.0/' + layer + '@EPSG%3A'+projection_epsg_no+'@pbf/{z}/{x}/{-y}.pbf'

I get the following error:

  GET http://192.168.10.15/geoserver/gwc/service/tms/1.0.0/test:poly@EPSG%3A4326@pbf/11/472/1228.pbf 404 (Not Found)

How do I pull in vector data from my geoserver using Openlayer's VectorTiles?


texts: Adding new rule by using QgsRuleBasedRenderer to existing one in PyQGIS?

texts:
I want to add new rules for a layer. However when I add them previous one(s) are deleted. Here is my code:

layer = iface.activeLayer()
symbol = QgsSymbol.defaultSymbol(layer.geometryType())
renderer = QgsRuleBasedRenderer(symbol)

# get the root rule
root_rule = renderer.rootRule()

# create clone of the default rule
rule = root_rule.children()[0].clone()

# set the label, expression and color
rule.setLabel("Akarsu Kutlesi")
rule.setFilterExpression(""""Subbasin" IN (2,4,6)""")
rule.symbol().setColor(QColor("orange"))

# append the rule to the list of rules
root_rule.appendtChild(rule)
# apply the renderer
layer.setRenderer(renderer)

Here is existing rules of layer:


texts: dissolve polygons of GeoJSON with gdal/ogr

texts:
I have a GeoJSON with only polygons in it. I have to dissolve (merge) the polygons with the same "DN"-value. I found out that you can do it with shapefiles like this:

ogr2ogr outputfile.shp inputfile.shp -dialect sqlite -sql “SELECT dissolvefield,ST_Union(geometry) as geometry FROM inputfile GROUP BY dissolvefield”

but since I dont know shp nor SQL I don't know how to apply that to GeoJSON and "DN"-value. I hope someone can help me out here.

{
    "type": "FeatureCollection",

    "features": [
    { "type": "Feature", "properties": { "DN": "0" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 1425, 0 ], [ 1425, 1 ], [ 1592, 0 ], [ 1425, 0 ] ] ] } },
    { "type": "Feature", "properties": { "DN": "1" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 4397, 88 ], [ 4397, 89 ], [ 4398, 88 ], [ 4397, 88 ] ] ] } },
    { "type": "Feature", "properties": { "DN": "2" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 369, 93 ], [ 369, 101 ], [ 371, 93 ], [ 369, 93 ] ] ] } },
    { "type": "Feature", "properties": { "DN": "2" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 368, 106 ], [ 368, 107 ], [ 369, 106 ], [ 368, 106 ] ] ] } }
    ]
}

EDIT:
I think I had a wrong Idea what "dissolve" is. So I want to explain what I want more detailed: I only want to MERGE two or more polygons with the SAME DN-Value AND which are neighbours. So that Polygons which were splitted due to tiles are unified again.


texts: How to convert small overlaping polygon features to inner ring in QGIS

texts:
I have a polygon vector layer with overlapping features (like below) and I would like to convert the features that are fully contained inside the bigger polygons (ie feature 2 and 4 below) into inner rings. I have tried saving the inner ring features as a separate polygon layer (they are all identified by attribute value) and using these as a difference layer in QGIS to no avail (I get a geometry error and result is only those features that were not intersecting). There a number of similar looking questions on here but I have not found any where the polygon vector layer contains distinct features to be made into inner rings.

Is there any way to just set those features as rings? I know I could digitise the rings but I need to do this for many many features and a more general programmatic solution would be much appreciated...


texts: Creating Random Points that respect minimal distance across features / strata using ArcGIS for Desktop?

texts:
I have a shapefile with multiple polygons in which I want to distribute random points. Each polygon should contain exactly one point. The minimal distance between the random points should be 200m.

Now the ArcGIS tool "Create Random Points" has an option "Minimum allowed distance". This distance is only respected for random points within a polygon / stratum (see 1). How can I distribute random points where the minimal distance is resptected across polygons / strata?

This question was already asked as How to generate random points with minimum distances between points in multiple polygons. However, in my view, the question was not sufficently answered.

Found in the ArcGIS 10.2 Help:

  This value will determine the minimum allowed distance between random
  points within each input feature. [..] Random points may be within the
  minimum allowed distance if they were generated inside or along
  different constraining feature parts.


texts:
I'm trying to present a mapserver with OL3, it doesnt work well. The documentation of the mapserver is non-existent and the GetCapabilities to the mapserver responds with a very short (not so good or usefull) answer.

I have managed to present the map, but it looks very weird, quite pixly. A picture says more then thousand words, so I'll just present to you how it looks.

http://imgur.com/UUVJEVK
http://imgur.com/lkSk3H4
http://imgur.com/mlNbHF0

And I even did a JS fiddle: http://jsfiddle.net/qpz7bp2g/

Any ideas?


texts: Getting extreme values for NDWI calculation for Landsat 8 Imagery

texts:
I am currently in a remote sensing course at my university and am working on a change detection project. I am having some trouble creating an ndwi raster from landsat 8 imagery. I am using arcgis 10.5.

I have downloaded the imagery from USGS earth explorer and am running the raster calculator tool to get my output. The formula I am using is:

float(band 3 - band 5)/float(band 3 + band 5)

My output contains low value of -20.4 and high values of 31.3 which definitely seems off to me. I ran a similar calculation on landsat 7 imagery from 2001 and got a result that ranges from -1 to 1. I am very new to the remote sensing field and am stumped on why my results for the landsat 8 imagery are so extreme.

Can anyone provide insight into why this is happening and what I can do to correct my calculations?


texts: Data copied from Excel to ArcMap layout view not exporting?

texts:
I prefer to use excel to make legends as I find it quicker than using ArcMap. Occasionally in some projects when I copy these in they appear fine in ArcMap, but some aspects don't export.

ArcMap layout view:

PDF export:


texts: MapServer: How to find all high and low spots in a raster?

texts:
I'm currently tinkering with MapServer to see if it could work as a backend for the remodel of my hobby stormchaser website.
Is it possible (in MapServer) to find all high and low spots in a raster grid and put a label on them. Obviously, I'm trying to make an atmospheric pressure map with "H" and "L" markers.
I know how to find the highest and the lowest value in a grid with GDAL but I basically want to find all points which are highest/lowest among all their neighbors.
A solution inside of MapServer would be preferred but I wouldn't be opposed to some preprocessing script.


texts: Get KML region info from address and kmz

texts:
I know nothing about GIS so I could be using the wrong words here.

I've got a KMZ file which shows the congressional districts of the US. I want to be able to enter a street address, and then get back which congressional district I'm in.

Can anyone give me some pointers on how to get started with this?


texts: How to enable the map unit function in "Sphere_Miller_Cylindrical" CRS in QGIS?

texts:
I am making  world maps, for varying purposes, and one of the requirements is to label each country with her name.
Normally I used the map map unit function in order to adjust the labels size to the map scale.
This time I was asked to use the "Sphere_Miller_Cylindrical" CRS, not the regular "WGS_84" CRS, so I change the project properties.

However, once I made this changed the label were gone.

What I need to do in order of making the label reappear?


texts: I would like to make the CartoDB heatmap use weighted values from a column in the dataset. I tried to find the option in the wizards and the docs but I couldn't find anything. Is there a way to do that?


texts: File stubbornly refuses to allow me to edit in QGIS

texts:
I opened up a shapefile ZIP in QGIS. Now, I cannot "toggle editing mode" to on, nor access any of the editing features, anymore! What could I be doing wrong?

Please, and also the one, here-for, this:


texts: Exporting symbols rather than points in .dxf

texts:
I have a map I've made in QGIS which is points from a .csv file. These have then be turned into different symbols based on what they represent and then the sizes of said symbols based on values in the .csv.

I need to export this as a .dxf to use in AutoCAD so it basically looks the same as it does in QGIS. However, when I do so it just exports as single dots (points). Is there a way I can export them as the symbols as I can see in QGIS?


texts: Testing layer for Has Selection using fidSet not working as expected

texts:
I am using ArGIS Desktop Standard 10.4.1 with Python 2.7.11, and am a S elf-taught Python scripter.

I've researched this extensively but there doesn't appear to be anything newer than 2014.

I am trying to test if a feature layer has a selection or not. I'm doing this via my own python script in ArcToolbox.

Essentially what is happening is that the length of len(arcpy.Describe(lyr).fidSet.split(";")) when no features are selected is returning 1. Prior posts on this subject indicate that it should return 0. Below is some very simple code:

lyr = arcpy.GetParameterAsText(0)    # parameter is a feature layer

arcpy.AddMessage(len(arcpy.Describe(lyr).fidSet.split(";")))
arcpy.AddMessage(arcpy.Describe(lyr).fidSet.split(";"))

if len(arcpy.Describe(lyr).fidSet.split(";")) > 0:
    arcpy.AddMessage("Has Selection")
else:
    arcpy.AddMessage("No Selection")

If zero feature are selected this returns

1
[u'']
Has Selection

If one feature is selected it returns

1
[u'5']
Has Selection

If more than one feature is selected it returns

9
[u'1', u' 2', u' 3', u' 4', u' 5', u' 6', u' 7', u' 8', u' 9']
Has Selection

I've tried testing the returned string with is None and also putting not in front of the length statement, but those still return Has Selection when no features are selected.


texts: Displaying intersecting roads and creating multiple labels for long roads in QGIS?

texts:
I am new to GIS and I am using QGIS 1.7.3 on Windows.

Can any one tell me how to display a map from line shapefile with these two characteristics:

Roads are displayed with borders at the sides and intersect at crossroads
Name label for long roads are duplicated at specified intervals

Here is the link to an image that illustrates my idea:

(The map was taken from underdark's blog)

This is the best I can achieve by drawing lines with two overlapping simple line, one slightly bigger than the other:


texts: QGIS 3.14 hangs when I use Data Source Manager on Windows 10

texts:
I am using 3.14.1, just installed. QGIS launches, and I can open a project. If I click the Data Source Manager button, it will appear to open the frame of a window, but then it hangs. It also displays "Getting Tiles" in the bottom bar of the main window.
Eventually the application crashes.
This was working in 3.12, and my upgrade seems to have broken it.
I have found someone else with this problem in this forum, but they are using QGIS 3.10.
Any ideas of a solution? Or what other information I should collect before I create a ticket?


texts: Tile error: Earth Engine capacity exceeded (exporting image) - sentinel2

texts:
I am new to earth engine. This is my second project, I followed eart engine tutorial on mygeoblog.com, when I try to export the image, it says ” Tile error: Earth Engine capacity exceeded”. I have tried to change the scale, but it still shows like that. Can you help me with that? I really need this map to do my landcover map.
This is my code: https://code.earthengine.google.com/4e22f5381df78c66d796d6d48b4c5a0e
This is my studyArea asset: users/dmdewi/Sumut3
Can you help me exporting this?
Thankyou!


texts: R gdal_grid interpolation results in null values for elevation

texts:
I am trying to grid point cloud data. I created a csv with coordinates x,y and elevation z-values.

The first lines of the .csv file look like this:

"X", "Y", "Z"
666578.24179169,5747856.3361029,-0.24
666578.251458881,5747856.34642544,-0.24
666578.250803538,5747856.36641516,-0.24
666578.250803538,5747856.36641516,-0.25
666578.260143058,5747856.38673256,-0.25

I run the following code in R:

gdal_setInstallation()
valid_install <- !is.null(getOption("gdalUtils_gdalPath"))

temporary_dir <- "directory"
tempfname_base <- file.path(temporary_dir,"dem")
tempfname_csv <- paste(tempfname_base,".csv",sep="")

# Now make a matching VRT file
tempfname_vrt <- paste(tempfname_base,".vrt",sep="")
vrt_header <- c(
  '<OGRVRTDataSource>',
  '\t<OGRVRTLayer name="dem">',
  '\t<SrcDataSource>dem.csv</SrcDataSource>',
  '\t<GeometryType>wkbPoint</GeometryType>',
  '\t<GeometryField encoding="PointFromColumns" x="X" y="Y" z="Z"/>',
  '\t</OGRVRTLayer>',
  '\t</OGRVRTDataSource>'
)
vrt_filecon <- file(tempfname_vrt,"w")
writeLines(vrt_header,con=vrt_filecon)
close(vrt_filecon)

tempfname_tif <- paste(tempfname_base,".tiff",sep="")

# Now run gdal_grid:
gdal_grid(src_datasource=tempfname_vrt,
                    dst_filename=tempfname_tif, a="invdist:power=2.0:smoothing=1.0",
                    txe=c(666360.2, 666580.4),tye=c(5747818, 5747878),outsize=c(400,400),
                    of="GTiff",ot="Float64",l="dem",output_Raster=TRUE)

The code does not give any error, but results in only null-values (elevation).


texts: QGIS project file saved to PostGIS – issue with users' privileges

texts:
I have multiple QGIS projects with PostGIS layers saved to different schemas (new feature in QGIS 3.x). I created various users with different privileges, but for the sake of this problem, let's assume that that there are 2 PostGIS users: one with all privileges and the other with read-only abilities.

When I connect to a PostGIS database using read-only credentials and load layers alone into QGIS, everything's fine, the user can only read and cannot edit layers in any way.

When I load the project, all layers and symbology are correctly loaded with it, but the read_only user appears to have all the privileges associated with the user that created the project. Even on the layers list when I move the cursor over the layer, it points to a different user.

At first I suspected QGIS might have somehow saved privileges of the admin user and use them to bypass PostGIS. I checked on other PC, on which I had never logged to any of the databases, and the problem persists.

Has anyone encountered such problem and perhaps have a solution? I wouldn't want to stop using that feature (saving QGIS projects to PostGIS), it's very convenient.


texts:
I'm using Server and Desktop 10.3. I'm trying to publish a geoprocessing tool as a service.This tool is a Python script that extracts raster data, outputs this as csv, and creates a chart in Excel using a VB script. This tool runs perfectly in Desktop. I've registered and validated all of the folders and data referenced in the script on the server. When I try to publsh, I get the error "

  {(u'Data: C:\Website_Test\Python_Scripts\Test_Scripts\Test_Charts\Build_raster_Charts.vbs, used by Script ChartTest cannot be copied to the server', 178): []}

This VB script is contained in of of the folders that has been registered with the server. According to the ESRI help page, the data store entry should cover this, yet it still won't work. Has anyone else had a similar issue?


texts: Error with Script Publishing Service from MXD with UNC Paths

texts:
I have script that attempts to publish a service which consists of a single tiff.

I am using boilerplate ESRI code for this one:

# initialize
imagery_date = 'some date'
folder_name = 'name of folder'
summary = 'summary info_{}'.format(imagery_date)
service_name = 'name of service_{}'.format(imagery_date)
workspace = r'\\xxxx\xxxx\xxxx\Imagery'
ags_con = r'xxxx\xxxx.xxxx.xxxx_6443 (admin)'
service = "{}/{}/{}.MapServer".format(ags_con, folder_name, service_name)
map_doc = arcpy.mapping.MapDocument(os.path.join(workspace, '{}.mxd'.format(service_name)))
sddraft = os.path.join(workspace, '{}.sddraft'.format(service_name))
sd = os.path.join(workspace, '{}.sd'.format(service_name))

analysis = arcpy.mapping.CreateMapSDDraft(map_doc, sddraft, service_name, 'ARCGIS_SERVER',ags_con, False, folder_name, summary, tags)

if analysis['errors'] == {}:
    # Execute StageService
else:
    # if the sddraft analysis contained errors, display them
    print("ERROR!!")
    print(analysis['errors'])

The script fails with this error:

ExecuteError: ERROR 001272: Analyzer errors were encountered (codes = 3).

This error is supposed to occur when the layers data source is not accessible, and ESRI suggests repairing the broken link.

When I open the MXD I'm publishing from, there are no broken links. When I analyze the service, there are also no errors, and I can publish from the MXD.

The tiff in my MXD is referenced via a UNC path, and I do not have relative paths stored.

Has anyone dealt with this issue before?


texts: Does OpenStreetMap have property boundaries?

texts:
I want to move from Google Maps to OpenStreetMap. I was planning on using MapBox to create the tiles if that makes a difference. OpenStreetMap seems to lack property outlines, which Google Maps has. I was wondering if this information is in OpenStreetMap but hidden, or if I would have to make my own custom layer.

A few screenshots to help show what I'm talking about. First Google Maps, then MapBox below.

So I would like to have the light grey lines and property numbers in the OpenStreetMap version. The map is for New Zealand. Thanks in advance for any help.

Update:
Turns out the data is freely available here:
https://data.linz.govt.nz/layers/category/property-ownership-boundaries/


texts: Improving fuzzy PNG image used as company logo in ArcMap?

texts: Why is my company Logo in the title block of my ArcGIS map blurry?  The logo was inserted into the title block as a PNG image and appears blurry when I export the map to PDF.


texts:
I'm trying to use CartoDB API to retrieve a dataset from my account to a heatmap. I want to filter the data using SQL. The problem is the following: the sql query does not always work when using the on done callback. However, it does work when I set it from the console using the map variable (where I stored the layer).

cartodb.createLayer(map, layerUrl)
.addTo(map)
.on('done', function(layer) {
    map.layer = layer;
    layer.getSubLayer(0).set(subLayerOptions); // sql is inside here
    map.invalidateSize(true);
}).on('error', function() {
    console.error("Cannot create layer");
});

EDIT

I wrote another cartodb code in the example, the actual sentence which sets the sql is this one: layer.setSQL(subLayerOptions.sql); (as it is a heatmap and the torque layer seems not to have the getSubLayer function).

To clarify my question, the result I'm obtaining is the same as if I don't set the query in the done callback (as I assume it is using the sql I have in the map I created from CartoDB SaaS). It is ignoring the query passed through javascript. I remark that if I set it later it works fine, it just needs some time to load the first time so I assume the problem is around that but I haven't find it yet.

EDIT 2

Sometimes the fix is just to zoom in and out. In this cases the data seems to be loaded but the interface is not refreshing.


texts: Seeking appropriate statistic to measure spatial autocorrelation of points with binary values

texts:
I am trying to determine the level of spatial autocorrelation in a point dataset.  The attribute I'm interested in is binary (presence/absence of a species), for which Moran's I is not appropriate.  On the other hand, Joint Count statistics, which are typically recommended for binary or categorical data, apparently are not appropriate for point data.
What is an appropriate statistic to measure global and/or local spatial autocorrelation of points when the attribute of interest is binary?


texts: Openlayers2 proxy.php cURL settings for GeoServer

texts:
We have to use a proxy to send a request to GeoServer. I created a PHP file like this.

OpenLayers.ProxyHost = "proxy.php?url=";

is Openlayers settings.

localhost/myapp/proxy.php

is my proxy adress. Proxy content is like this.

<?php
$headers = array(
    "Content-type: application/xml",
    "Content-type: text/xml",
);

 $ch = curl_init($_GET['url']);
 curl_setopt($ch, CURLOPT_HEADER, 1);
 curl_setopt($ch, CURLOPT_HTTPGET, 1);
 curl_setopt($ch, CURLOPT_FOLLOWLOCATION, 1);
 curl_setopt($ch, CURLOPT_RETURNTRANSFER, 1);
 curl_setopt($ch, CURLOPT_HTTPHEADER, $headers);
 $output = curl_exec($ch);
 curl_close($ch);
 echo $output;
?>

I am using this proxy like:

http://localhost/myapp/proxy.php?url=http://somedomain.com

Which works, such that a web page is opened.

But when I try a GeoServer WMS get query, it gives a NullPointerException error. For example:

http://myip/geoserver/app/wms is returning admin page, but

http://localhost/myapp/proxy.php?url=http://myip/geoserver/app/wms is returning:

HTTP/1.1 200 OK Server: Apache-Coyote/1.1 Content-Type: application/xml Transfer-Encoding: chunked Date: Wed, 03 Sep 2014 14:18:11 GMT java.lang.NullPointerException null

http://myip/geoserver/app/wms?request=GetCapabilities is returning an XML page, but

http://localhost/myapp/proxy.php?url=http://myip/geoserver/app/wms?request=GetCapabilities

is returning plain text result.


texts: Tool in QGIS to Get Area Weighted Average

texts:
I am looking for a tool that can give me an area-weighted average.

I have a shapefile of parcels and another shapefile with Curve Numbers. I want to get Area Weighted CN for each parcel. If I do it manually, I would first intersect parcels with the Curve Number file. Recalculated Area for each feature in the intersected shapefile. Take the attributes to Excel and then use the Pivot table or Vlookup to find the sum of areas and sum of (CN*Area of each feature) for a particular parcel and then divide them to get the Area Weighted Curve Number.

Is there a tool to do it automatically in QGIS?


texts: Line pattern fill in legend symbology- QGIS 2.18

texts:
I work with "line pattern fill" in polygon layer with distance value: 5

and the legend is:

My problem is when i choose value 20 in the distance option- the green line disappear from the symbology legend.

How can i choose distance value 20 and still the green line will be presented in the legend? i know one solution by draw a green line in the legend by myself, but i looking for a smart automatic solution.


texts:
Having tested out the functionality of OpenGeo Suite Community, I have installed the Linux version on my Centos 5 Server. However unlike in windows the dashboard is greyed out for Post GIS. I believe this is because the local web interface limits access due to security? where as in windows you can bypass this by loading the desktop client. I don't think you can load the local GUI without using a browser on Centos?

I read I need to run PGadmin separately, however I cannot find it in Linux when it should have installed along with the OpenGeo Suite Community?

In short how can I access the PGadmin interface? or re-install it? I found the Enterprise GUI installers do not work similarly to other posters on here.


texts:
I've done an intersection with the ITopologicalOperator6 interface, and I now have an IGeometry with the result.

How to draw this in a graphics layer on the map? I want to see what was intersected (possibly colored).

According to this topic: CustomLayer Draw Com Exception

it seems that I could create a new layer from the IMap interface.. is this possible? Is there any other more convenient way to do this?

How would you draw an IGeometry item to your map?


texts: Extracting columns from shapefile into new raster using ArcGIS Desktop?

texts:
Using ArcGIS 10.4

I have a line (shapefile) with lots of information. In 8 columns there are groups of species (Group 1 to Group8), indicating if this is the dominant species group in the area or not. Then for every row there is one group that has a 1 (dominant) and all the other 7 have a 0. In rare cases non of the species groups is dominant (all are 0).

Now, I need ONE raster file with this information. My idea was a raster with the numbers 1 to 8 indicating the dominant species group. But I cant get it to work. First I used feature to raster to extract each of the 8 columns. Then I tried to merge them into one raster with the mosaic to new raster tool. But of course this only gives me a file with 0 and 1 again.

How could I get this done?


texts: Publish Kernel Density Layer Which is Transparent

texts:
I am using ArcGIS Pro 1.4 to publish web layers which are the result of a kernel density operation.  I am using mining data from MRDS for China.  The raster layers were created and published to AGOL here: http://tiles.arcgis.com/tiles/nAQVsl698blIdU3d/arcgis/rest/services/Major_Mineral_Extraction_in_China_WTL1/MapServer.

When they are activated in either application (Pro or AGOL) all but one of the layers is not fully transparent in the areas of the published extent which do not contain the output cells.  To see this, turn on all of the layers at once.  Only the layer KernelD_MRDS_TI is transparent in those areas.  This prevents proper visualization of the combined layers when they are all active.  As a result, the slight grey tint accumulates and blurs out the entire basemap.  I created all of the layers with the same tool and with the exact parameters with the exception of the population field used.  Here is the vector point data used in the kernel density tool: http://services1.arcgis.com/nAQVsl698blIdU3d/arcgis/rest/services/Major_Mineral_Extraction_in_China_WFL1/FeatureServer.

The clip here shows the parameters used for the kernel density tool:

I would like all of them to be as fully transparent in those regions.


texts: Animating quantitative values in Carto Torque Cat?

texts:
In Carto (CartoDB), is it possible to change the category classifications in “Torque Cat” to be qualitative instead of categorical, as you would have in the “Choropleth” map? Both the values and the legend representation. Basically I want to animate the representation we get from the choropleth wizard.

I’m working on a simple animated map of points over time for drone locations and I want the animation to change the point and trail colors according to a field with values for the drone roll or pitch at each location. I've attempted to port over the CSS from the choropleth map to the torque cat map but not getting the results I want.

updated to add my modified code as applied from answer below:

/** torque visualization */

Map {
-torque-frame-count:512;
-torque-animation-duration:30;
-torque-time-attribute:"serialdate";
-torque-aggregation-function:"sum(pitch)";
-torque-resolution:2;
-torque-data-aggregation:linear;
}

#seacharger_ca_hi_10kmjoin{
  comp-op: source-over;
  marker-fill-opacity: 0.9;
  marker-line-color: #FFF;
  marker-line-width: 0;
  marker-line-opacity: 1;
  marker-type: ellipse;
  marker-width: 6;
  marker-fill: #0F3B82;
}
#seacharger_ca_hi_10kmjoin [ value <= 12] {
   marker-fill: #41b6c4;
}
#seacharger_ca_hi_10kmjoin [ value <= 6] {
   marker-fill: #7fcdbb;
}
#seacharger_ca_hi_10kmjoin [ value <= 0] {
   marker-fill: #ffffcc;
}
#seacharger_ca_hi_10kmjoin [ value <= -7] {
   marker-fill: #FFCC00;
}
#seacharger_ca_hi_10kmjoin [ value <= -14] {
   marker-fill: #FF9900;
}
#seacharger_ca_hi_10kmjoin [ value <= -21] {
   marker-fill: #FF6600;
}


texts: Using is_in with area or relation ID in Overpass API

texts:
I am trying to find a way to determine the admin_level hierarchy of an area (e.g. island) in Overpass. I understand that I can easily find all sorts for things within the area, but I need to find out where a particular area "is_in".
Edit:
I can easily find admin_levels within the area using this query (thanks TomazicM for the link):
[out:json];
area[name="Isola d'Elba"]->.island;
nwr[admin_level](area.island);
out;

But what I need is the opposite - all admin_levels in which Elba is in.
Is this possible at all?


texts: arcmap layer export as shapefile, in the dbf file the time value is missing

texts:
In arcmap, I export a layer as shapefile. After that I find from the dbf file, the editor tracking date-time columns (created_date, last_edited_date) value only contain the date value, The time value is not available.

Is it normal? How can I get the time value in the dbf file also?

I am using arcmap 10.2.2.


texts: MapInfo: How to get the population contained within 2 geographical tables

texts:
I using 3 different tables which are also geographical objects:

population data; points on the map representing inhabitants (1 or more)
region data; the map is divided in regions
coverage data; it contains several coverage areas overlapping or not

What I am searching for is: What is the percentage of people covered by a certain coverage area in a certain region?

I tried several things but nothing worked, for example with SQL select I tried:

Select Columns: Regions.name, CoverageArea.name, Sum(Population.inhabitants)

from Tables: Regions, CoverageArea, Population

where Condition: Regions.obj Contains Population.obj And CoverageArea.obj Contains Population.obj

Group by Columns: 1, 2

Order by Columns: 1, 2

into Table Named: Selection

"And I get Incorrect tables are joined. Invalid join condition in WHERE clause."

Hope I was clear enough :s

** edit **

Few tests I did. I can't understand why the thrid SQL query doesn't work...


texts:
I am trying to develop a web app like this. That is, I want to show a map with points and when the users clicks on a point, it should display some information about it.
Can I do that with TinyOWS?
This app requires WMS, WFS or both services?


texts: Point to Polygon nearest distance DS_distance is not using geography index & knn  or  does not give result in order

texts:
The below query gives result in wrong order.

SELECT *, ST_Distance_Sphere(ST_SetSRID(ST_MakePoint(-75.66,45.33), 4326 ) , geometry ) / 1609.34 AS distance
  FROM myTable
 ORDER BY geometry <#> ST_SetSRID(ST_MakePoint(-75.66,45.33), 4326 ) LIMIT 1

Index on the geometry & geography exists.
I have also tried st_distance in order by, results are perfect but It does not use geometry or geography index.
I have tried st_DWithin but it is also taking time. for a shorter distance, Dwithin is working fine, but in this case, search distance is at least with 5000 miles.

PostGreSQL version: 9.4
No of rows: almost a million.

how knn fits in above scenario?

I have upgraded version from 9.4 to 9.5. It's working Updated working query is:

SELECT *, ST_Distance_Sphere(ST_SetSRID(ST_MakePoint(-75.66,45.33), 4326 ) , geometry ) / 1609.34 AS distance
FROM myTable
ORDER BY geometry <-> ST_SetSRID(ST_MakePoint(-75.66,45.33), 4326 )::geography LIMIT 1


texts: Using iterate function to extract time series in Google Earth Engine?

texts:
The code below extracts average rainfall for each US county for one given day and send the table result to drive. Now, how can I use iterate function to extract a time series for one year for example?

 // load precip data (mm, daily )
var precipCollection = ee.ImageCollection('OREGONSTATE/PRISM/AN81d')
                    .select('ppt')   // select  precip band only
                    .filterDate('2017-01-01'); // e.g. one day
print(precipCollection, 'precipCollection');

//reduce the image collection
var dailyPrecip = precipCollection.reduce(ee.Reducer.sum());
print(dailyPrecip);

// visualize annual precipitation
var precipPal = ['white','blue','red','pink'] // store palette as variable
Map.addLayer(dailyPrecip, {min: 1, max: 100, palette: precipPal}, 'precip');

////////////////////////////////////

// load regions: counties from a public fusion table, removing non-conus states
// by using a custom filter
var nonCONUS = [2,15,60,66,69,72,78] // state FIPS codes that we don't want
var counties = ee.FeatureCollection('ft:1ZMnPbFshUI3qbk9XE0H7t1N5CjsEGyl8lZfWfVn4')
        .filter(ee.Filter.inList('STATEFP',nonCONUS).not());
print(counties, 'counties');

// visualize
Map.addLayer(counties,{},'counties');

// get mean precipitation values by county polygon
var countyPrecip = dailyPrecip.reduceRegions({
  collection: counties,
  reducer: ee.Reducer.mean(),
  scale: 4000 // the resolution of the prism dataset
});
print(countyPrecip);

// select features
var out = countyPrecip.select(['.*']);

// add a new column for date
out = out.map(function(feature){
return feature.set('date','2017-01-01');
});

// Table to Drive Export Example

Export.table.toDrive({
  collection: out,
  description: 'precip_by_county',
  folder: 'GEE',
  fileFormat: 'CSV'
});


texts: Generating graph edge list from polyline and point shapefiles?

texts:
I have a polyline shapefile of rivers/stream segments, and a point shapefile of intersection points (nodes) between the segments (edges). Each stream segment starts and finishes at an intersection point. I wish to generate an edge list file (.csv or similar) which contains the start and end node IDs for each stream segment ID, in a format similar to "Edge ID, From ID, To ID" (see example image below).

Preferably I would like a solution to this within ArcGIS 10. So far I've experimented with a iterative approach (iterating through each edge segment and selecting nodes within a +-5m radius), but feel that this might not be optimal for large stream networks (100,000+ intersections).

I have uploaded a basic sample dataset (one polyline and point shapefile matching the image above) here.


texts: Making Google Maps tiles of EPSG:4258 maps

texts:
I have a map image, EPSG:4258, which as I found out is within bounds: -10.67W 71.05N 31.55E 34.5S.

gdal_translate -of Gtiff -a_ullr -10.67 71.05 31.55 34.5 -a_srs EPSG:4258 file.png file.tif
gdalwarp -s_srs "EPSG:4258" -t_srs "EPSG:3857" file.tif file-converted.tif
gdal2tiles.py --profile=mercator -z 1-8 file-converted.tif output

It should cover the area of the UK, so it's a bit too north.
At another website here's what they ended up with:

Should I use some conversion like cs2cs for the a_ullr attribute?


texts: Removing black line from bar diagram in QGIS 3.14

texts: I am trying to remove or move the black "line" which you can see on the some countries but I can't do anything.


texts:
I am trying to merge some shapefiles with different schemas in GeoTools. I'm aware of ogr2ogr, but for now I would like to figure out this problem with GeoTools.

The steps I'm taking are:

Create a new Merged shapefile with the schema from the first shapefile in a list
Update the schema of the Merged shapefile file with the subsequent shapefiles in the list of shapefiles
Add features to the merged shapefile
Write the shapefile

The problem is, I can't add features from shapefiles that contain different schemas. OK, that seems straight forward enough.

But when I try to get the schema from a subsequent shapefile and update the Merged shapefile schema, GeoTools throws an exception:

"Schema modification not supported"

Why won't GeoTools let me update the schema? Is there something more than just simply calling

mergedData.updateSchema(shapefileData.getSchema())

public static void init(ArrayList<String> shpArray, String sessionName)
{
    Config config = new Config();

    String sessionData = config.getProperty("pathTo_fade_servlets_data") + "sessions/" + sessionName;
    String newUrl = sessionData + "/tempShp.shp";

    init(newUrl, shpArray.get(0));

    for(int k = 1; k < shpArray.size(); k++)
    {
        mergeShapefile(shpArray.get(k));
    }

    saveShapefile();
}

@SuppressWarnings("rawtypes")
public static void init(String newURL, String copyURL)
{
    try
    {

        File newFile = new File(newURL);

        Map<String, Serializable> newFileParams = new HashMap<String, Serializable>();
        newFileParams.put("url", newFile.toURI().toURL());

        ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();

        mergedData = (ShapefileDataStore) dataStoreFactory.createNewDataStore(newFileParams);

        //

        File copyFile = new File(copyURL);

        Map<String, Serializable> copyFileParams = new HashMap<String, Serializable>();
        copyFileParams.put("url", copyFile.toURI().toURL());

        ShapefileDataStore shapefileData = (ShapefileDataStore) dataStoreFactory.createNewDataStore(copyFileParams);

        mergedData.createSchema(shapefileData.getSchema());

        mergedFeatureStore = (FeatureStore)mergedData.getFeatureSource(shapefileData.getTypeNames()[0]);

        mergedTransaction = mergedFeatureStore.getTransaction();

    }
    catch (MalformedURLException e)
    {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
    catch (IOException e)
    {
        // TODO Auto-generated catch block
        e.printStackTrace();
    }
}

@SuppressWarnings({ "rawtypes", "unchecked" })
private synchronized static void mergeShapefile(String url)
{
    if (url == null) return;

    try
    {
        File file = new File(url);

        Map<String, Serializable> fileParams = new HashMap<String, Serializable>();
        fileParams.put("url", file.toURI().toURL());

        ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();

        //

        ShapefileDataStore shapefileData = (ShapefileDataStore) dataStoreFactory.createNewDataStore(fileParams);

        FeatureCollection features = shapefileData.getFeatureSource().getFeatures();
        //error occurs here:
        mergedData.updateSchema(shapefileData.getSchema());

        mergedFeatureStore.addFeatures(features);

    }
    catch (Exception e)
    {
        e.printStackTrace();
    }
}

 private static void saveShapefile()
 {
        if (mergedTransaction == null) return;

        try
        {

            mergedTransaction.commit();

            mergedTransaction.close();

        }
        catch (Exception e)
        {
            e.printStackTrace();
        }
    }


texts: Binding DOM Elements to Graphics Instead of SimpleMarkerSymbol()

texts:
Can you let me know if it is possible to bind our custom DOM into the geometry point instead of ArcGIS API SimpleMarkerSymbol()?

I tried to do this like below but it is not doing the job!

var point = new Point(-126.687482, 54.793577);
var pointSymbol = '

';
var pointAttributes = {city: "Albuquerque",  state: "New Mexico"};
var pointInfoTemplate = new InfoTemplate("Albuquerque");
var pointGraphic = new Graphic(point, pointSymbol, pointAttributes).setInfoTemplate(pointInfoTemplate);
    thePoints.push(pointGraphic);
        console.log(thePoints)
map.on('load', function () {
        for (i = 0; i < thePoints.length; ++i) {
        map.graphics.add(thePoints[i]);
    }
});


texts: Seeking detailed topographical map of Istanbul?

texts:
I need to acquire a very detailed topographic map of Istanbul, Turkey without any markers. Do you know such sources who provides these kind of maps?  I tried to Google it, but the results were far from satisfactory. I found world topography from NASA, some low resolution images that does not worth much. There are a lot for US but unfortunately I failed to find one for Istanbul.

Actually since I don't the details of GIS, I don't know the technical terms to explain what I am looking for. But here is what I'll use the map for: I am going to print 3D version of Istanbul using some 3D printer, and display visualizations on top of it using projector hooked on the ceiling. In this case, I assume that I need either a 3D point cloud, or a 2D height map where I can deduce 3D points after some processing.


texts: Accessing GPS data with JavaScript in real time

texts: Is it possible to access GPS data in real time on a smart phone with JavaScript (without an android or ios application I mean)?


texts: How to accumulate in the image collection in Google Earth Engine

texts:
Here is an image collection [img1,img2,img3,img4,img5]. I want to get List or image collection [img1, img1+img2,img1+img2+img3, img1+img2+img3+img4, img1+img2+img3+img4+img5].

I tried to complete the following code https://code.earthengine.google.com/20763e4a37f2341558c2b7162f02eb83. But the result image collection b did not contain 'system:time_start'. I don't know how to modify my algorithm. Or you have a better way to get the accumulated image collection with 'system:time_start'.

var img1 = ee.Image(1).set('system:time_start','2001-01-01');
var img2 = ee.Image(2).set('system:time_start','2001-01-02');
var img3 = ee.Image(3).set('system:time_start','2001-01-03');
var img4 = ee.Image(4).set('system:time_start','2001-01-04');
var img5 = ee.Image(5).set('system:time_start','2001-01-05');
var image=ee.ImageCollection([img1,img2,img3,img4,img5]);
/*image = image.map(function(image) {
  return image.add(1);
});*/
var list = image.toList(image.size());

var t = ee.List(image.map(function(tt){
  return tt.get('system:time_start');
}))

var newList = list.slice(1);
newList = newList.iterate(function(data, list) {
  data = ee.Image(data);
  list = ee.List(list);
  list = list.add(data.add(ee.Image(list.get(list.size().subtract(1)))));
  return list;
}, ee.List([list.get(0)]));
var b = ee.ImageCollection(ee.List(newList).slice(1));

print(b)


texts: Automatically give colours to shapefiles in different layers

texts:
I'm in the final stage of my floodmap and I want to colour the River basin districts based on the Total floods in the area. When there are 0 floods it has to be white, 1-2 floods -> green, 3-5 -> yellow and go on to red. For example all the river basins from Austria has their own attribute tabel (not all in one table). as you can see below.

I already made a colourstyle which i would like to use. It is working now if I would do it manually, but I would like to have it automatically because in the future the floods will count on. For example if there is one flood, it could be 3 or 4 in 2020. Now the colour is green because of the one flood, but with 3 or 4 it has to be yellow.

Later on I am going to link my Excel to arcmap (have no idea how, but that's a later problem), so when I add a flood in excel it will automatically change from attributes and colours in arcmap.

Does somebody know a solution for this problem?
(the give colours based on attributes to different layers problem)


texts: Making raster layer 1 meter high using ArcGIS Pro

texts:
I have 30 raster layers of the subsurface. Each representing the subsurface 1 meter below one another. I would like to have a 3D view of this like voxels. Each layer must have an artificial height of 1 meter.
There is no need to make calculations with the volume. Just the picture "must look nice".
How do I do this using ArcGIS Pro?


texts:
I am trying to display my bing map raster tiles in QGIS. I am using an HTML file like this, except with my own raster tiles, as described in WMS Tile Layer Example.

My own raster tiles don't have a projection, they have quad keys.

Is there a way to add this to QGIS with Python?

I am guessing it is something to do with this, but I cannot quite figure it out. Got through some documentation PyQGIS Cookbook | Raster Layers but with not much success.

My data is currently stored on S3, I have tried this:

from osgeo import gdal

xml = """<GDAL_WMS>
    <Service name="VirtualEarth">
           <ServerUrl>
                  https://s3.amazonaws.com/bucketname/key/t_${quadkey}.png
           </ServerUrl>
    </Service>
    <MaxConnections>4</MaxConnections>
    <Cache/>
</GDAL_WMS>"""

vfn = "/vsimem/bing.xml"
gdal.FileFromMemBuffer(vfn, xml)
rasterLyr = QgsRasterLayer(vfn, "sm")
rasterLyr.isValid()
QgsProject.instance().addMapLayers([rasterLyr])


texts:
How do I decipher these cryptic COMExceptions that ArcObjects are throwing?

I got this exception System.Runtime.InteropServices.COMException: HResult: 0x8004151C and thats it, no more info no nothing. How do I go about debugging?

I managed to decipher the exception after hours of googling, talking to colleagues, searching ESRI's forums and resources and at last gis.stackexchange.com.

I found this resource ArcObjects error codes and some other resource I cant manage to find again, where I was supposed to use the calculator in scientific mode, and do 10 steps of stuff get a 10 digit error code.

While doing so I couldn't help but get annoyed and wonder if this is the way its meant to be, or if there is any solution, code or anything else out there that could ease the job.

Or do we developer really have to wait for an exception to happen, and then do a manually lookup on the internet and maybe, just maybe find the cause for the exception after a long period of time?

So thats why I'm asking here, hence the question in the top.


texts:
I have a points feature collection and have added the latitude and longitude to each point's properties as well. Now I want to do a processing on each point and want to use the map function. The processing function will extract some Landsat 8 time series and do an analysis, and the analysis subfunction should be fed with the parameters in a dictionary format and will return a dictionary.

def getModel(pointFeature):
    args = pointFeature.toDictionary()
    lon = args.get('longitude')
    lat = args.get('latitude')
    point = ee.Geometry.Point([lon, lat])
    # Create image collection
    l8_collection = ee.ImageCollection(
        'LANDSAT/LC08/C01/T1_SR').filterBounds(point).filterDate('2015-01-01', '2018-01-01')
    info = l8_collection.getRegion(point, 30).getInfo()
    header = info[0]
    data = np.array(info[1:])
    results = analysis_function(header, data)
    return ee.Feature(None, results)

Results = points.map(getModel)

But I get the below error and it traces to the line info = l8_collection.getRegion(point, 30).getInfo():

ee.ee_exception.EEException: Failed to decode JSON.
Error: Field 'value' of object '{"type":"ArgumentRef","value":null}' is missing or null.
Object: {"type":"ArgumentRef","value":null}.

I think it is all about client/server variable passing in GEE, but I couldn't manage to solve it. I tried to add getInfo() after calculations of Lat and Lon or casting the Lat and Lon to ee.Number() but the problem persists.
Any help?


texts: Creating raster from selection of footprints in Mosaic Dataset?

texts:
I am working with thousands of DEM tiles, so I have put them into a Mosaic Dataset. (Please give me feedback if this is not the best way to handle this amount of data).

I was given DEM tiles that extend beyond the area of interest, so I did a Select by Location in order to choose the footprints that fall within the area of a particular polygon.

I have tried right clicking the "Image" title in the Mosaic Dataset and choosing Export data (with Selected Footprints). This returned a raster with all of the rasters in the Dataset.

EDIT
I have tried Exporting the data from the Mosaic Dataset twice. Once using a Mask and again using Extents (all in Geoprocessing Environments). For both I used a polygon I created. Yet the output is still the full raster.

Background geoprocessing is turned off, so I can't figure out what is going wrong so far.

How do I go about creating a raster from of the area of interest?


texts: Cannot display Oracle raster images using QGIS

texts:
I am dealing with QGIS 2.2 and the plug-in Oracle spatial georaster.

Apparently the connection to the database can be made, as the raster images are listed once you press 'Connect', but when try to display them, either I see a black box or anything at all.

In theory my QGIS 2.2.0 has the GDAL version 1.10.1. I downloaded it using the OSGEO4w package for Windows.

Is there any way of testing where the problem might be?