texts:
I am considering trying to make an API so that I can make some spatial datasets available to colleagues for analysis.

Part of my work has been to analyze and prepare data which can then be used for further analysis by others. The work (while currently at a smaller scale and less sophisticated) is similar to walkscore but does involve some enormous datasets. There are increasing restrictions on how I can share the original data, but my derivative work is shareable. I have been thinking about about how best to share the results of my analysis (outside of passing on large datasets) and thought that an API would be one approach. What sort of things, should I be thinking about when constructing an API? Are there design specifications that I can follow?

My vision sounds a little more grandiose than it currently is, but I think it would be a useful framework to consider early on in this work.


texts: Get the inverse of default heat map gradient in Google Maps Javascript API

texts: In Google Maps Javascript API for heat maps the default gradient displays red color in areas of higher density and changes to red as intensity decreases. In my application I need to show areas of higher intensity in green and lower intensity in red. I tried printing the default gradient values heatmap.get('gradient') only to get 'undefined' value. How can I get the default value of Google Maps API heat map gadient.


texts: Unable to install OSM Editor on ArcGIS 10.2.0?

texts:
Esri have provided installers of the ArcGIS Editor for OpenStreetMap for 10.2.1 and 10.2.2, but not for 10.2.

How can I make it?


texts: Can I convert one big .bmp file (size is 19GB (high resolution city map)) into  one or several .tiff format files without dividing this .bmp file into smaller tiles?


texts: Override typographic ligatures in QGIS labels?

texts:
For one of my label layers I am using a font that contains typographic ligatures for certain letter pairs. If I increase the letter spacing in Layer Properties, the ligature is retained as a single character, so the spacing between the pair of letters does not increase (see image for example—note the double T ligature).

Is there a Python expression, or another way to disable or override these ligatures, so that the individual characters are displayed instead?


texts: Unexpected behaviour of RasterBand converted to Array and back

texts:
I have two raster images, mask and data. The first one is a single-band tif with integer values that refer to classes. The second one, data, is a 3-band (RGB) tif. In the raster containing the classes, I assigned the value 100 to pixels that could not be identified.
Now I tried to set the R = G = B = 0 in data for every pixel where mask has the value 100.
To this end, I wrote the following code:
(For easier understanding of variable names: a = array, b = band, m = mask)
import numpy as np
from osgeo import gdal
mask = "C:/Users/Manuel/Nextcloud/dat/omk/03_2021/ortho_B3_4_mask.tif"
data = "C:/Users/Manuel/Nextcloud/dat/omk/03_2021/ortho_B3_4_CROP.tif"
NoDataVal = 100

from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
def blacken(data, mask, NoDataVal, outFile = "/vsimem/Blackened.tif"):
    dat = gdal.Open(data, GA_ReadOnly)
    msk = gdal.Open(mask, GA_ReadOnly)
    bm = msk.GetRasterBand(1)
    am = BandReadAsArray(bm)
    null = np.where(am == NoDataVal)
    am = None
    bm = None
    msk = None
    drv = gdal.GetDriverByName("GTiff")
    new = drv.Create(outFile, dat.RasterXSize, dat.RasterYSize,
                     dat.RasterCount, dat.GetRasterBand(1).DataType)
    CopyDatasetInfo(dat, new)
    for b in range(dat.RasterCount):
        band = dat.GetRasterBand(b+1)
        a = BandReadAsArray(band)
        a[null] = 0
        new.GetRasterBand(b+1).WriteArray(a)
        a = None
    new.FlushCache()
    dat = None
    return new

However, the values that are set to 0 are not the ones I asked for. It looks like the removed pattern is upside-down. Maybe this is due to the images location in the southern hemisphere? Or some issue in my code?
I attached an image that compares the 3 layers: The original data, the mask layer where the value 100 is the large gray matrix surrounding the patches, and the output image, where the black area is supposed to resemble the gray matrix of the mask.
Obviously, the patches don't fit while everything else seems to line up. Looking at the larger pattern, the black area appears to be vertically mirrored. The images were previously cropped to the same size and they are in the same reference system. I thought I treated them in the same way, so I wonder why the blackened pixels don't match the pixels where mask = 100. Perhaps there is an annoyingly simple mistake somewhere...


texts: Use coordinates from Geotools geometry-object to send request to Webserver

texts:
I have a polygon as a Geotools Geometry Object.

No I want to use the coordinates of the polygon to send a request to a webserver.
The webserver demands the coordinates to be in the following format (all coordinates separated by commas):

http://....&coordinates=21.205285624773502,62.75338836781611,21.241110826971651,62.75063258303164,21.231190001747548,....

Is there a name for that kind of format and has Geotools a method to convert it directly to it?

The only way I see right now would be to convert it to a String, which gives me a WKT-format and then using String methods to get it to the right format. I thought there might be a better way.


texts: Specifying LAS dataset filters in ArcGIS ModelBuilder?

texts:
I'm trying to create a tool in the ArcMap ModelBuilder to extract a canopy height raster from LiDAR data. I need to create rasters for ground elevation and first return. I've got my LAS dataset, and the "Filter" option on the LAS dataset toolbar allows me to select whether to display/use the ground points, first return points, etc.

I want to be able to do this in ModelBuilder - to be able to automatically create rasters for ground and first return. But I can't figure out how to specify these filters in the LAS Dataset to Raster tool or in the properties of the input LAS datasets themselves.

Does anyone know how to do this?


texts: Labeling a joined (one to many) layer with the number of feature joined to it in ArcGIS Online

texts:
Using ArcGIS Online, I performed a one to many join with a buildings layer (one) and a table of internal cameras (many). This produced a layer showing buildings and allows me to scroll though all the cameras associated with a particular building.
Now, I would like to label the buildings based on the number of cameras joined to a building. For example, a building with 8 cameras would be labeled 8.
I know I will have to populate a field with the number of cameras in each building, but I am not sure how to do this. Any recommendations?


texts: v.what.rast this tool does not work in qgis 3.2?

texts:
I have tried to use the v.what.rast tool in Qgis 3.2 to extract elevation values from an image of the alos-palsar sensor. And it does not generate the column with the elevation data. I create the elevation column previously to be filled by the tool.
In terms of the fields, it is mentioned in the GRASS v.what.rast help page that the field type must be numeric:

Integer, Float, Double

So you will need to create a numeric field, not string, in order to store the results.

I follow the steps of the tool in Qgis and I always get the table of the new empty shapefile.

Why is this problem.

I prepare the tool. and I run the algorithm

This is the result I get ...
What is the problem?

WARNING: Missing column name

Update vector attributes ...

0..WARNING: There is no record for category 3 in the table

14..WARNING: There is no record for category 2 in the table

28..WARNING: There is no record for category 4 in the table

42..WARNING: There is no record for category 1 in the table

57..WARNING: There is no record for category 5 in the table

71..WARNING: There is no record for category 0 in the table

85..WARNING: There is no record for category 6 in the table

100

v.what.rast complete. 0 updated records.


texts: How to change the raster cell values where they overlap with the cells of another raster?

texts:
I have a raster that I'm trying to edit by changing the value of select cells (those that overlap with the cells of another raster).

For some context, what I'm trying to do is create a Distance to Water layer (Euclidean distance) where the interior of water bodies is "NoData". This has all been derived from a watercourse polyline layer and a waterbody polygon layer. I'm using ArcGIS 10.3.1.

Below is a screenshot the two rasters I'm working with. The green is the waterbody layer, and I'm hoping to change the value of the cells in the Distance To layer that overlap to NoValue (or something else which I can later reclassify as NoValue).


texts: ModelBuilder appends .shp to output file after running

texts:
I'm running model that processes line features into a raster, and I'm using a raster to point tool to further sieve the result.

When the model hasn't run yet, these are the tool parameters:

As you can see, the result is as a GDB feature class, and everything work perfectly.

Yet after running the tool once, the output now has a .shp appended to it, and since it's written in a gdb, it obviously throws an error.

Is there a way to override this or is it a bug?

This is the only tool that does so in a model with more than 150 tools.


texts: How to prevent creating blocks when using DWG writer to draw lines or polygons in FME desktop

texts: I had a MapInfo TAB reader and a DWG writer in my workspace. I used line part of DWGStyler transformer to draw lines and polygons in my DWG output file, when I ran my workspace the output schema was correct but FME created lots of blocks for most of the lines and polygon. Is there any way to prevent creating blocks?


texts: When is data created in datums other than WGS84?

texts:
I'm trying to understand why we need datums other than WGS84. Since GPS coordinates are relative to the WGS84 ellipsoid what is the advantage of having other datums like for example NAD83. I understand that other datums are locally a closer approximation to the earths topological surface - but when is this advantageous in real life?

I feel the answer may be in this term from the wikipedia entry for Geodetic Datum, but still unclear how you would create location references to NAD83 without a GPS.

"A reference datum is a known and constant surface which is used to describe the location of unknown points on the Earth. Since reference datums can have different radii and different center points, a specific point on the Earth can have substantially different coordinates depending on the datum used to make the measurement. There are hundreds of locally developed reference datums around the world, usually referenced to some convenient local reference point. Contemporary datums, based on increasingly accurate measurements of the shape of the Earth, are intended to cover larger areas. The most common reference Datums in use in North America are NAD27, NAD83, and WGS 84.

The North American Datum of 1927 (NAD 27) is "the horizontal control datum for the United States that was defined by a location and azimuth on the Clarke spheroid of 1866, with origin at (the survey station) Meades Ranch (Kansas)." ... The geoidal height at Meades Ranch was assumed to be zero, as sufficient gravity data was not available, and this was needed to relate surface measurements to the datum. "Geodetic positions on the North American Datum of 1927 were derived from the (coordinates of and an azimuth at Meades Ranch) through a readjustment of the triangulation of the entire network in which Laplace azimuths were introduced, and the Bowie method was used." (http://www.ngs.noaa.gov/faq.shtml#WhatDatum ) NAD27 is a local referencing system covering North America.

The North American Datum of 1983 (NAD 83) is "The horizontal control datum for the United States, Canada, Mexico, and Central America, based on a geocentric origin and the Geodetic Reference System 1980 (GRS80). "This datum, designated as NAD 83 ...is based on the adjustment of 250,000 points including 600 satellite Doppler stations which constrain the system to a geocentric origin." NAD83 may be considered a local referencing system."


texts: Reason for "L.Control.Draw is not a constructor" error?

texts:
Below is my code:

<!DOCTYPE html>

        <title>Simple Leaflet Map</title>
        <meta charset="utf-8" />
        <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7/leaflet.css"/>
        <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/leaflet.draw/0.2.3/leaflet.draw.css" />

I don't know why I'm getting this error and also L is not defined error is coming while adding leaflet.draw.js.


texts: Removing the first lines of an ascii file as a python script tool

texts:
I want to delete the header of an ASCII file with a python script in ArcGIS. The header is always six lines long.
After running the script that needs this ASCII file, I need to reinsert the header. The header is always the same text.
At the end it will be part of a Python-script-tool for the ArcGIS Modelbuilder.
Can anyone help me?


texts: Export/import from local PostGIS db to cartodb?

texts: I have a local PostGIS database with a polygon table and a nonspatial table. I would like to upload them into my CartoDB account. What is the best way to do this? I could export the polygon layer as a shp and the nonspatial table as a csv then upload them like that, is that the most straightforward way? Or is there a way to skip that and copy directly from postgis to postgis?


texts: Alternative to MapTiler for tiling a GeoTIFF

texts:
I'm currently using MapTiler Pro console version, called from C# code as process with the following parameters:

  -geodetic -nodata 0 0 0 -f png32 -o result/ -tms -resampling bilinear -srs "+proj=longlat +datum=WGS84 +no_defs" -zoom 8 12 source.tif

Is there any alternatives to MapTiler for using from code? I tried vips's dzsave, but output differs from needed and it lacks some functions, like tiling exact levels.

I also use gdal's c# bindings, but it lacks gdal2tiles methods, which, from what I've found, is the only alternative.


texts: How does Resampling raster work at corner cells?

texts:
This may seem quite a simple question for which I couldn't find a satisfactory answer. I have been trying to understand how resampling of raster work with ArcGIS 10.3, in both cases of refining and coarsening, especially at corner and edge cells.

For descriptive purposes, I have chosen an input raster with 4 cells of size 1x1m which I resample to a cell size of 0.25x0.25m using different resampling techniques.

With bilinear resampling technique, I got this result:

What I couldn't understand here is why the four cells at all four corners have the original values. If it is taking distance weighted average of four input cells, the values would be different. Or ArcGIS is adding three extra cells (of same value as the corner cell) at the corners beyond the edges of the input raster such that four closest cells of the four corner cells on out raster always have the same value of original corner cell? Something like this:

One more thing to mention is when I used cubic method, I got the same result as with bilinear method. In this case, the input raster has only four cells; how does cubic method calculates a distance averaged value from 16 cells?


texts: What's the QGIS Legend Format Syntax to have thousands separated by commas?

texts: I have numbers in the millions and would like to know what string to enter in the Legend Format box to have thousands separated by commas ex: 2,000,000


texts: How to search text string and select item in QGIS?

texts:
I have a table attribute with complete roads name.
I have to search items by digiting a part of the road name and I need to colour all the items selected.

eg: ROAD NAME: / road 01 / road 02 / road 11 / road 12

I need to colour all items that contain "1" (so all the items except the second one).
I would like to set the string operator is this way:
ROAD NAME ILIKE ' * 1 * ' but it does not work

edit I work in advanced search in table attribute


texts: Examples of leaflet interactive mapping applications

texts:
I would like to develop a simple mapping application. One where users can digitize lines, polygons and points and attribute the features. The attributes will then be stored in a database for later assessment.

Has anyone seen such an app using the leaflet library?


texts:
I am having problems installing Anaconda with ArcGIS 10.1. My problem is specifically with pandas. ArcGIS python is my default python version. I installed python 0.16.2 in this version, which installed numpy 1.9.2. Following some posts, I created a .pth file (with ZZ as preffix, in order to read first the default version modules) and put it in C:\Python27\ArcGIS10.1\Lib\site-packages, which points to the Anaconda packages (C:\Anaconda2\Lib\site-packages). I can import most of the Anaconda python modules, but I have a problem with pandas. When I import pandas, I get the following error:

RuntimeError: module compiled against API version a but this version of numpy is 9

This only happens in a Windows command window, in ArcGIS python window imports the 0.16.2 version. The Anaconda pandas and numpy versions are 0.17.1 and 1.10.1. When I import the numpy version in a command window I get 1.9.2 (ArcGIS installed numpy version). The only problem is the error above with pandas, but after the error I check the version and says 1.9.2. I have read about creating environments in Anaconda, but I need the pandas 0.16.2 version because of some new funcionalities not present in older versions. I tried uninstalling the ArcGIS pandas version to use the Anaconda version, but I got the following error:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Anaconda2\Lib\site-packages\pandas\__init__.py", line 7, in <module>
    from pandas import hashtable, tslib, lib
  File "pandas\src\numpy.pxd", line 157, in init pandas.hashtable     (pandas\hashtable.c:38262)
ValueError: numpy.dtype has the wrong size, try recompiling

Does anyone how can I solve this problem?


texts:
I'm currently trying to get a couple of layers uploaded to a geoserver and I'm having trouble getting the projection recognized. The result of the curl command is that the two layers are both missing the Native SRS, which should be EPSG: 3857.

I'm able to fix by manually going to the layer, choosing a Declared SRS, and computing the Lat/Lon Bounding Box from the native bounds. This allows me access to the layer via WMS ,which is what I want.

But this is not a good solution. I have a lot of layers that I need to upload, so I don't want to manually fix the missing Native SRS problem for each layer. So how can I fix the shapefiles so that geoserver can find the Native SRS on its own?

I got the correct .prj file from here: http://hub.qgis.org/issues/7368
So both the .prj and the .qpj files are set to this for the two shapefiles:

PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],AUTHORITY["EPSG","3857"]]


texts: Geocoding Result layer does not show up in ArcMap?

texts:
http://i.imgur.com/XF5iQN4.png

Everything was going fine until I geocoded the addresses. The addresses were all matched, but the Geocoding Result layer wouldn't show itself on the map. There should be about 40 points that are missing.
Any ideas what I'm doing wrong?
I already checked for coordinate systems; they're the same.


texts: Deleting collinear vertices from polygon feature class using ArcPy

texts:
I have a shapefile with polygons and want to write the Python script to eliminate collinear vertices in a polygon. I found a solution for the line feature class, but not for polygon: Deleting collinear vertex / vertices using ArcMap.
import arcpy

#----------------------------------------------------------------------
def are_collinear(p1, p2, p3, tolerance=0.5):
    """return True if 3 points are collinear.
    tolerance value will decide whether lines are collinear; may need
    to adjust it based on the XY tolerance value used for feature class"""
    x1, y1 = p1[0], p1[1]
    x2, y2 = p2[0], p2[1]
    x3, y3 = p3[0], p3[1]
    res = x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)
    if -tolerance <= res <= tolerance:
        return True

#----------------------------------------------------------------------
def get_redundant_vertices(vertices):
    """get redundant vertices from a line shape vertices"""
    indexes_of_vertices_to_remove = []
    start_idx, middle_index, end_index = 0, 1, 2
    for i in range(len(vertices)):
        p1, p2, p3 = vertices[start_idx:end_index + 1]
        if are_collinear(p1, p2, p3):
            indexes_of_vertices_to_remove.append(middle_index)

        start_idx += 1
        middle_index += 1
        end_index += 1
        if end_index == len(vertices):
            break
    return indexes_of_vertices_to_remove

#----------------------------------------------------------------------
def clean_geometries(fc, densify_curves=False):
    """clean polyline features in the fc removing redundant vertices"""
    in_sr = arcpy.Describe(fc).spatialReference.factoryCode

    with arcpy.da.UpdateCursor(fc, ['SHAPE@', 'OID@']) as ucur:
        for row in ucur:
            print "OBJECTID", row[1]
            cleaned_parts = []
            shape = row[0]

            if 'curvePaths' in shape.JSON:
                if densify_curves:
                    shape = shape.densify('DISTANCE', 1, 1)
                else:
                    continue

            for part in range(shape.partCount):
                vertices = [(p.X, p.Y) for p in shape.getPart(part)]
                if len(vertices) < 3:  #polyline's part consists of 2 vertices
                    continue
                vertices_to_remove = get_redundant_vertices(vertices)
                vertices_to_keep = [
                    val for idx, val in enumerate(vertices)
                    if idx not in vertices_to_remove
                ]
                cleaned_part_as_array = arcpy.Array(
                    [arcpy.Point(*coords) for coords in vertices_to_keep])
                cleaned_parts.append(cleaned_part_as_array)

            if cleaned_parts:
                cleaned_shape = arcpy.Polyline(
                    arcpy.Array(cleaned_parts), in_sr)
                row[0] = cleaned_shape
                ucur.updateRow(row)

if __name__ == '__main__':
    fc = r'C:\GIS\Temp\ArcGISHomeFolder\Default.gdb\_SimpleRoads'
    #r'C:\GIS\Temp\ArcGISHomeFolder\empty.gdb\roads'
    clean_geometries(fc, densify_curves=False)

If I use this code for a polygon FC it will delete my polygons. I am new to Python and I don't know how to fix it in ArcPy.
Could you give me some tips or ideas how to achieve it?


texts: I am looking to use realcentroid where I'm assuming it puts all points within the polygons rather than the tool in QGIS Vector -> geom -> centroid tool, which in my result some bendier polygons have their centroids placed outside - the geometric centroid? The real centroid plugin doesnt install in QGIS v3 - do these get updated with releases?


texts: how to allow null values in a feature-layer

texts:
I'm using ArcDesktop 10.2. I have 30 x,y locations that have 20 fields of numeric data, stored in an Excel spreadsheet. Some cells have no data, and I want those cells to be represented as null values when I convert it from a spreadsheet to a feature. I keep getting 0's instead of < NULL > when I convert the Excel spreadsheet to a shapefile. The ArcGIS help website says that in version 10.2, they no longer support null values in shapefiles. When I try to change the field properties regarding "Allow Null Values", it will not let me change the setting from No to Yes.

It's critical for a statistical analysis later on that the null values remain and are not converted to 0's. Is there a way to create a feature class that allows null values? There are some true 0's in the dataset, so I cannot simply delete all the 0 values using the Editing Toolbar. I could compare the shapefile and the original spreadhsheet cell-by-cell and manually remove the 0's that are actually nulls, but I figure there must be a better way.


texts: SAGA tool polygon self-intersection: time to prosess increases exponentially with increasing buffer width

texts:
With QGIS 3.10 I try to find overlapping areas and the number of overlaps of buffers with the SAGA tool "Polygon self-intersection". The tool works fine with round buffers around points with a buffer-width of five kilometers. When I tried to find the overlapping area with ten kilometer buffers around the same points the tool takes forever while qgis does not respond anymore. To approximate the problem I widened the buffer from five to six to seven kilometers and checked the processing time in the report. This was 6, 9 and 90 seconds respectively. So the wider the processing time increases exponentially.
To solve the problem I tried:

to reduce the number of buffers, but there was no difference.
buffers in Shape-File vs buffers in geopackage. But also no difference.

Can anyone confirm this behavior?
Unfortunately I also need the number of overlaps per overlapping area. So the tool is kind of perfect. Is there another Open Source alternative? So far I only know of the tool "find overlaps" in ArcGIS Pro.


texts: Building a ol3 web application with dynamic content and user registration

texts:
My set-up:

Currently I’m using a PostGIS database, with a GeoServer service and displaying the data with OpenLayers 3.  I use WFS data to display the features in clusters and I’m able to edit them with WFS-T.  The PostGIS database is very large ~250.000 features but I load them with the WFSBBOXLoader strategy. With the ol3edit template from Boundless, this works very well.

My question:

What is the best way/set-up to allow only registered user to modify features with the ol3 WFS-T system, such that the user can register themselves?

Should I even use my current set-up or other API's?
I just need the direction/ideas I should go... :)

PS: I don't like drupal+extensions...


texts: Extracting coordinates from multipoint geometries?

texts:
I found this guide to get pixel value of raster under point: Getting pixel value of GDAL raster under OGR point without NumPy?

I've a problem here:

for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

If I print geom, I can see all the point's coordinates but if I print mx or my the result is 0 for every feat.

if I print geom:

MULTIPOINT (9.48556682927651 39.1368546561669)
MULTIPOINT (9.46173837787432 39.1298387934389)
MULTIPOINT (9.54058251881085 39.441437960101)
MULTIPOINT (9.08688113454087 39.5390542536278)
MULTIPOINT (9.00215692841357 39.0072908447031)
MULTIPOINT (9.17741005408364 39.2751570664614)

if i print mx or my:

0.0
0.0
0.0
0.0
0.0
0.0

I'm using Python 3 maybe this could be a problem


texts: Reversing Extract by Mask in ArcGIS for Desktop?

texts:
I have a DEM and I want to apply some tools to it (fill and others). In my DEM I've zones (drawn as mask) that need to be exlcluded from that processing.

I don't want to apply the processing to the whole zone and then put a mask on it.

How do I get the zone where I can apply my tools? My masks are the zones I don't want to process.

I thought about Clip (data management) but it does not work, I also tried to substract this mask from the DEM, without success. Same with "Erase tool", it does not accept my DEM.

I'm using ArcMap 10.2.2.


texts: Trying to clip raster in QGIS, "Cutline is not Valid"

texts:
I am using QGIS 2.14.5-Essen and I am trying to simply clip an image using a square polygon. I am encountering an error that I have never seen before. This is what I am attempting to do:

"clipper" is a shapefile containing a rectangular polygon feature over my desired study area.

The image I am trying to clip is a LANDSAT-8 Surface Reflectance image. When I load it into QGIS, it is already in spatial reference system

+proj=utm +zone=35 +datum=WGS84 +units=m +no_defs

and I verified that it is showing up in the correct area by loading in the OSM basemap and confirming that it is oriented properly.

My "clipper" shapefile is also in

+proj=utm +zone=35 +datum=WGS84 +units=m +no_defs

However, when I try to run the operation,which corresponds to the gdal command

gdalwarp -q -cutline D:/african_towns_identifying/data/original_image/clipper.shp -tr 30.0 30.0 -of GTiff D:/african_towns_identifying/data/original_image/LC81720712016212LGN00_sr_band2.tif D:/african_towns_identifying/data/clipped_image/band2.tif

I get this error:

I don't understand why it just won't clip normally, because I have done this before without any problems. I don't understand the error at all and cannot find anything helpful with a Google search

Trying to clip the image manually using the "Extent" in the dialog box and dragging a box will execute the tool successfully, but I get an image ranging from NaN to NaN


texts: In CartoDB, is it possible to associate metadata with dataset columns?

texts:
I have some data in google spreadsheets. I've been moving the data to cartodb.com because I want to use CartoDB for mapping. As well as mapping, I want to be able to access the raw data and present it using C3 charts.

Moving away from a generic spreadsheet to an SQL-backed system, I've hit a limitation. I want to be able to show friendly names and tool-tip descriptions for each metric displayed in a chart, but this is data that applies to a whole column.

I can do this in google spreadsheets by reserving specific rows for metadata, but importing this in to cartodb does not work for columns that contain non-string values.

I thought about having a secondary dataset in CartoDB that contains metadata for a primary dataset. Of course, this means making two requests to get all the data for a chart but I assume this would not be prohibitive.

Are there any other workarounds to accomplish this?


texts: Obtaining extents of each data-driven page using ArcPy?

texts:
I have a set of index polygons (rectangles) that vary in size and orientation on which data-driven pages are based. I am looking for a method that would give me extents Xmin, Xmax, Ymin, & Ymax for each page.

Is there a built-in way to do this, perhaps with the ArcPy mapping module or can I look at the centroid of each index polygon, do some math based on page size, and then have the desired extents?

EDIT:
I was playing around with this and I think I've started to locate a solution. Using the dataframe's extent property I am able to get the information I need. Now, all I have to do is iterate over the data driven pages and get the extent of each.

mxd = arcpy.mapping.MapDocument(Map_Path)
DF = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
for page in pages:
  print DF.extent


texts: Accessing First Member Coordinates of geoJSON

texts:
I apologize if this is trivial or not possible, I am new to GIS. So I have a geoJSON file with thousands of members, each with a point coordinate. Is there a way for me to access the first coordinate (and the last, but I will deal with that later) of the geoJSON file?

Extraneous info: I am working on an html page using leaflet and have converted a CSV file into a geoJSON file, but I don't want the initial set view to be hard coded. I want it to be set, ideally to the center of the points, but being set to the first point will also work.


texts: Pass geometry as input to custom PostGIS function

texts:
I am desigining a custom PostGIS function and I am at a loss on how to pass a geometry column from a table as an input to another column.

I am writing a function to return certain columns from a table, and overwrite their cost_s and reverse_cost_s values if they intersect an input geometry.

This is what I have so far:

CREATE OR REPLACE FUNCTION replace_cost(
    in_geom GEOMETRY
)

    RETURNS TABLE(
        class_id INTEGER,
        name VARCHAR(255),
        cost_s DOUBLE PRECISION,
        reverse_cost_s DOUBLE PRECISION
    ) AS

    $BODY$
        SELECT
            a.class_id AS class_id,
            a.name AS name,
            a.cost_s AS cost_s,
            a.reverse_cost_s AS reverse_cost_s
        FROM
            ways AS a
        WHERE
            /*a.class_id != test_input*/
            NOT ST_Intersects(in_geom, a.the_geom)

        UNION ALL

        SELECT
            a.class_id AS class_id,
            a.name AS name,
            'INFINITY' AS cost_s,
            'INFINITY' AS reverse_cost_s
        FROM
            ways AS a
        WHERE
            /*a.class_id = test_input*/
            ST_Intersects(in_geom, a.the_geom)

        ORDER BY
            cost_s

    $BODY$

    LANGUAGE SQL;

SELECT * FROM replace_cost(overwrite_poly.geom);

But I get the error back:

ERROR:  missing FROM-clause entry for table "overwrite_poly"
LINE 57: SELECT replace_cost(overwrite_poly.geom);

I think I kind of understand, because I don't have overwrite_poly in any kind of FROM clause, but I don't know how or where to fit it in.  I also can't find much in the way of writing functions that take a PostGIS geom as input.


texts: ArcGIS Flex application - toggling service layers with AS3?

texts:
I'm working on an application that has to access data from a server. The flex application uses Esri's World_Imagery basemap and is able to access my data from the ArcGIS for Server, running locally on the same machine.

There are multiple layers to the data that are be displayed, and I would like to toggle them from ActionScript. What is the best way to do this? I've taken a look at the documentation, but it's unclear to me how to go about doing this without the preset widgets that Esri provides.

Here's my code:

<?xml version="1.0" encoding="utf-8"?>
<s:Application xmlns:fx="http://ns.adobe.com/mxml/2009" xmlns:s="library://ns.adobe.com/flex/spark" xmlns:esri="http://www.esri.com/2008/ags" pageTitle="Example - dynamic maps on top of projected map tiles">

    <fx:Script>
        <![CDATA[
            var b:Boolean = false;
            if (b) {
                //Show Pennwell_PowerPlants
            } else {
                //Hide Pennwell_PowerPlants
            }
        ]]>
    </fx:Script>

    <esri:Map level="3" wrapAround180="true">
        <esri:center>
            <esri:WebMercatorMapPoint lon="0" lat="0" />
        </esri:center>
        <esri:ArcGISTiledMapServiceLayer url="http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer" />
        <esri:ArcGISDynamicMapServiceLayer url="http://localhost:6080/arcgis/rest/services/Pennwell_Substations/MapServer" />
        <esri:ArcGISDynamicMapServiceLayer url="http://localhost:6080/arcgis/rest/services/Pennwell_PowerPlants/MapServer" />
    </esri:Map>
</s:Application>


texts: Listening current tool Changed in ArcObjects10.1?

texts:
I tried listening current tool changed to show current tool's message and caption in my main form. but I find esri interface about nothing.

How to listen for the current tool being changed in ArcObjects 10.1 or are there any methods I can use to get current tool?


texts: Update 1 random value in a group of values

texts:
I have an attribute table with a date field and a random field. The date field has three different times for each day. I need to input a '1' in the random field for one of those three times in each day.

I created a field that removed the time stamps and just included the day in order to group each day. This is indicated as 'Date' in the code below.

This code sample inputs '1' in he Random column for each of the six time of a random day instead of just one for each day.

import arcpy, random
from arcpy import env
env.workspace = r"H:\Working_Data.gdb"

for fc in arcpy.ListFeatureClasses('*'):

  object_ids = [r[0] for r in arcpy.da.SearchCursor(fc, ['Date'], )]
  sample_size = len(object_ids) / 6
  random_ids = random.sample(object_ids,sample_size)

  with arcpy.da.UpdateCursor(fc, ['Date','Random']) as Cursor:
    for row in Cursor:
        if row[0] in random_ids:
            row[1] = 1
        else:
            row[1] = 0
        Cursor.updateRow(row)
        print str(row) + 'updating'

I have tried multiple solutions from similar questions but can't get anything to work.

Data:

Date_Time                Date
2003-12-12 1:00:54 AM    2003-12-12
2003-12-12 11:00:54 AM   2003-12-12
2003-12-12 4:02:17 PM    2003-12-12
2003-12-13 1:01:30 AM    2003-12-13
2003-12-13 11:00:24 AM   2003-12-13
2003-12-13 4:00:12 PM    2003-12-13
2003-12-14 1:00:06 AM    2003-12-14
2003-12-14 11:03:11 AM   2003-12-14
2003-12-14 4:00:32 PM    2003-12-14
2003-12-15 1:00:52 AM    2003-12-15
2003-12-15 11:01:54 AM   2003-12-15
2003-12-15 4:01:22 PM    2003-12-15

The output I'm looking for would be something like this where each day has one '1' placed randomly in the random column.

Date_Time                Date          Random
2003-12-12 1:00:54 AM    2003-12-12    1
2003-12-12 11:00:54 AM   2003-12-12    0
2003-12-12 4:02:17 PM    2003-12-12    0
2003-12-13 1:01:30 AM    2003-12-13    0
2003-12-13 11:00:24 AM   2003-12-13    0
2003-12-13 4:00:12 PM    2003-12-13    1
2003-12-14 1:00:06 AM    2003-12-14    1
2003-12-14 11:03:11 AM   2003-12-14    0
2003-12-14 4:00:32 PM    2003-12-14    0
2003-12-15 1:00:52 AM    2003-12-15    0
2003-12-15 11:01:54 AM   2003-12-15    1
2003-12-15 4:01:22 PM    2003-12-15    0


texts: Images overlay in 3D surface plot matplotlib (python)

texts:
It is possible to overlay a geotiff, png or jpg to a surface_plot generated in matplotlib?

I want to display for example a smooth hillshade, geology or orthoimages. Another question... how to smooth the surface result?

This is the code I used to generate the plot:

#open terrain data using gdal
dtm = gdal.Open('AR.asc')
gt =dtm.GetGeoTransform()
banda1 = dtm.GetRasterBand(1)
xres = gt[1]
yres = gt[5]

# The bounding box of a shapely geometry
#used to filter a portion of digital terrain model

bb = geom.bounds
minx,miny,maxx,maxy = bb[0], bb[1], bb[2], bb[3]

inv_gt = gdal.InvGeoTransform(gt)
_x0, _y0 = gdal.ApplyGeoTransform(inv_gt, minx, miny)
_x1, _y1 = gdal.ApplyGeoTransform(inv_gt, maxx, maxy)
x0, y0 = min(_x0, _x1), min(_y0, _y1)
x1, y1 = max(_x0, _x1), max(_y0, _y1)
zi = banda1.ReadAsArray(int(x0), int(y0), int(x1-x0), int(y1-y0))

fig = plt.figure()
fig.set_size_inches(18.5, 10.5)
ax = fig.add_subplot(111, projection='3d')
plt.hold(True)
X = np.arange(gt[0], gt[0] + zi.shape[1]*xres, xres)
Y = np.arange(gt[3], gt[3] + zi.shape[0]*yres, yres)
X, Y = np.meshgrid(X, Y)
ax.plot_surface(X, Y, zi, alpha=.9, linewidth=0,cmap=plt.cm.gist_earth,antialiased=False,
                shade=False,rstride=1, cstride=1)
plt.show()


texts: Why do base and overlay only match at scale 1:1000 and 1:2000?

texts:
We are using GEOSERVER 2.1.4 to  overlay WMS on GoogleMaps V3 with Openlayers .

The problem we have is at different zoom levels the scale of the WMS overlay is different from the scale of the base layer.

It is only at scale 1:1000 and 1:2000 that overlay and base layer coincide.

Any suggestions?

The part of the code:

        var map;
        var untiled;
        var tiled;
        var pureCoverage = false;
    OpenLayers.Lang.setCode('el');
        // pink tile avoidance
        OpenLayers.IMAGE_RELOAD_ATTEMPTS = 5;
        // make OL compute scale according to WMS spec
        OpenLayers.DOTS_PER_INCH = 25.4 / 0.28;

        function map_reg_init(){
            // if this is just a coverage or a group of them, disable a few items,
            // and default to jpeg format
            format = 'image/png';
            if(pureCoverage) {
                document.getElementById('filterType').disabled = true;
                document.getElementById('filter').disabled = true;
                document.getElementById('antialiasSelector').disabled = true;
                document.getElementById('updateFilterButton').disabled = true;
                document.getElementById('resetFilterButton').disabled = true;
                document.getElementById('jpeg').selected = true;
                format = "image/jpeg";
            }
    var projEGSA=new OpenLayers.Projection("EPSG:2100");
    var projLONLAT = new OpenLayers.Projection("EPSG:4326");
        var projGOOGLE=new OpenLayers.Projection("EPSG:900913");
            var bounds = new OpenLayers.Bounds(
                2327363.1529655424, 4805801.76429202,
                2328463.0753220413, 4806793.842097636
            );
    var scales = [500, 1000, 2000, 4000, 5000, 10000, 15000, 25000, 50000];
    var resolutions = [];
    for(var i = 0; i < scales.length; i++) {
        resolutions.push(OpenLayers.Util.getResolutionFromScale(scales[i], 'm'));
    }

            var options = {
                controls: [],
                maxExtent: bounds,
        scales: scales,
                maxResolution: 156543.0339,
            numZoomLevels    : 22,
                allOverlays      : true,
                projection: "EPSG:900913",
                units: 'm'
            };
    var mousePointerStyle = 'default';


texts:
The radius of gyration is a measure of the standard deviation of distances between points (and their center of mass), i.e., how frequently and how far points are moving.
A formal definition is:

I have a large table of localization data (x-y coordinates, not a geometry field - might be transformed if faster?). What is the best way to compute the radius of gyration with PostGreSQL / PostGIS?

For the center of mass, ST_Centroid(ST_Union(ST_MakePoint(x,y)))?
Then ST_Distance? How to aggregate all this? ST_Distance(MakePoint(x,y),ST_Centroid(ST_Union(ST_MakePoint(x,y)))) does not work (ERROR: column "x" must appear in the GROUP BY clause or be used in an aggregate function).

Sorry, maybe a stupid question, I'm not so familiar with PostGreSQL.


texts: Best practise for managing shapefiles/symbology in QGIS

texts:
One aspect of GIS that wasn't covered in my training was the best way to manage shapefiles within a project.

In my hypothetical example I have a shapefile with a number of attributes. I want to produce a map comparing 4 different combinations of various attributes.

So I want 4 map windows in my print composer, I then have three options for how to set up the layers (that I have thought of):

A) 1 layer, 1 shapefile. I then use the 'lock layers' and 'lock styles' for each  map window, adjusting the properties and styles of the shapefile for each desired attribute combination as I go. The problem is, there is no record of what rules or attributes I used if I come back to it in 6 months time.

B) 4 layers, 1 shapefile duplicated in each layer. I can set the propoerties and style for each layer to correspond with each submap but the shapefile itself still refers to the original .shp file.

C) 4 layers, 4 shapefiles. Duplicate the shapefile and load each separate shapefile into its own layer and set style and properties for each shapefile . This seems to cause performance issues for me.

So far, I tend to follow B, this lets me fine tune settings and styles with ease and also gives reasonable performance but I don't know whether this may give rise to other problems down the line. What is generally considered best practise in this situation for performance and stability?


texts: Slope analysis on Terrain data using WhiteBox GAT?

texts:
How do I find out the Slope of "best fit to the sectoral terrain" using WhiteBox GAT?

Here are the details about the area of calculation.

We have a reference point, from the reference point our sectoral area could be around 1000m with the sector range starts from 220 degree to 300 degree. I wanted to calculate the "best fit to the sectoral area between 220 degree to 300 degree with the length of 1000m".

What I did so far is, I calculated the slope of the whole terrain but I can't limit the slope area within 1000m with 220 deg to 300 degree. Again I don't know how to calculate the "best fit" of the slope for the sectoral plane.

This helps to finish this one.


texts: Building raster with multiple attributes based on 3 input rasters?

texts:
I am trying to make a new raster based on three input rasters which have the same cell size and the grids align perfectly.

I now want to construct a new raster, which for each cell, has three attributes, being the cell the values of the different input rasters. I've read about using 'combine', but this provides with an attribute table with a count of all cells with the same combination of values in the different input rasters (I hope this is clear). What I need is an attribute table in which all cells are individually represented.

Does anyone have any advice on this?


texts:
i tried leaflet search control. Searching works great, i can see list with my  points names, but when I want search it, it does not work, not move to search point. My JSON data are points.   `    var data = us_states;

var map = new L.Map('map', {zoom: 5, center: new L.latLng([37.8, 16]) });

map.addLayer(new L.TileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png')); //base layer

var featuresLayer = new L.GeoJSON(data, {
        style: function(feature) {
            return {color: feature.properties.color };
        },
        onEachFeature: function(feature, marker) {
            marker.bindPopup('
'+ feature.properties.OBJECTID +'
');
        }
    });

map.addLayer(featuresLayer);

var searchControl = new L.Control.Search({
    layer: featuresLayer,
    propertyName: 'OBJECTID',
    marker: false,
    moveToLocation: function(latlng, OBJECTID, map) {

        var zoom = map.getBoundsZoom(latlng.layer.getBounds());
        map.setView(latlng, zoom); // access the zoom
    }
});

searchControl.on('search:locationfound', function(e) {

    //console.log('search:locationfound', );

    //map.removeLayer(this._markerSearch)`


texts: Problem of projection with overlaying layers (GeoServer) on map (OSM) using Leaflet

texts:
I'm working on a web mapping application using GeoServer and Leaflet, the problem I have is when I overlay the WMS layer (loaded from GeoServer) on base map I do not have superposition.

This is my code:

var map = L.map('map', {
            layers: [OpenStreetMap], /// fond de base
            center: [30.412647, -9.579837],/// coordonnées
            zoom: 13    //// zoom par defaut
});

/////layers de base
var baseLayers = {
        "Open Street Map": OpenStreetMap,
        "World Imagery": WorldImagery
};

var dist_layer = L.tileLayer.wms("http://192.168.250.198:8082/geoserver/test_layer/wms", {
        layers: 'test_layer:Districts_Agadir',
        format: 'image/png',
        transparent: true,
        version: '1.1.0',
        attribution: "myattribution"
});

The result I get:

This is my problem exactly the layer does not match base map

My layer is a shapefile created using ArcGIS 10.2 with Projection:

Projected Coordinate System:    Sud_Maroc_Degree
Projection: Lambert_Conformal_Conic
False_Easting:  500000,00000000
False_Northing: 300000,00000000
Central_Meridian:   -5,40000000
Standard_Parallel_1:    29,70000000
Scale_Factor:   0,99961560
Latitude_Of_Origin: 29,70000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_Merchich_Degree
Datum:  D_Merchich
Prime Meridian:     Greenwich
Angular Unit:   Degree

This is how my layer looks in ArcGIS it matches the base map exactly, and that's what I'm looking for in my web mapping application.


texts: Caculating Longitude and Latitude for Moving Object

texts:
Possible Duplicate:
  Calculating Lat/Lng X miles from point

I don't know how to ask the question, so I'll explain my requirments -- I hope I'm in the right site :)

Let's assume that a person is standing somewhere on Earth and his coordinates are: x long, y lat.
Let's assume that if the person is facing north he is on 0°, south is 180°, north-west is 315° ... etc. (unless there is some standard that I'm not aware of)
Let's assume that the person can walk in a straight line from his position on the direction he's facing until he completes a full round around the Globe; i.e. return to his starting position.

My question is, what is the equation/set of equations to be used to calculate the person's coordinates (longitude and latitude) every 50, 100, 300 ... etc. kilometers until a full round is done?

Thanks in advance.


texts: Is there a established color convention for US Federal Agencies?

texts: Do US Federal Agencies have a convention for coloring the lands under the jurisdiction of each? For example, BLM lands are often yellow, USFS light green, NPS blue or dark green, etc. I've found this, but it seems to be for the internal use of a single agency, where I'm looking for public applications.


texts:
In QGIS 3.4, I would like to place random points inside each polygon using a field from an attribute in the polygon layer and tag each point with its polygon id.  Currently, I use Random Points inside Polygon.  The output is what I am looking for exempt that the function doesn't carry over the polygon ID.  Therefore, I have to do a Spatial Join back with the Polygon layer.  This step is taking a long time due to performing an linear search.

Is there a way to drop random points and track the polygon ID back to the point?  A plugin?


texts: Geodatabase feature to KMZ - display issue with stacked lines

texts: I am working with a line feature in ArcGIS. My end product is a KMZ for use in Google Earth. However, there are lines that are stacked on top of each other in the feature. When I convert my feature to a KMZ only one feature is clickable in Google Earth, meaning clicking on stacked lines only displays the attributes of the topmost line. Is there any sort of method to deal with stacked lines? I know that point features that are stacked expand when you click on them to display all stacked points, but this is not a feature with lines apparently.


texts: Extrapolating point values to raster, excluding NAs

texts:
I have a raster generated from drone imagery that has been classified into water (river channel) and NA (all areas outside the wetted channel). I have a point file running down the center of the river channel and in some cases along the side.

I would like to extrapolate the point values across all of the 1's, to create a new raster where the point values have been extended, but I do not want the extrapolation to extend across the NAs.

How might I do this in R efficiently?

In the photo provided, the colored points represent different types of channel units (pool, riffle, rapid, etc) and the black area represents the wetter channel or 1's.

Structure of my points file:

Classes ‘sf’ and 'data.frame':  5588 obs. of  9 variables:
 $ CU_Number : num  5 5 5 5 5 5 5 5 5 5 ...
 $ Fish_Site_: Factor w/ 1 level "UpperLemhi1_2019": 1 1 1 1 1 1 1 1 1 1 ...
 $ Site_ID   : Factor w/ 1 level "UpperLemhi1_2019": 1 1 1 1 1 1 1 1 1 1 ...
 $ StreamName: Factor w/ 2 levels "Lemhi River",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ CU_Type   : Factor w/ 9 levels "LSC pool","LSC Pool",..: 9 9 9 9 9 9 9 9 9 9 ...
 $ Seg_Number: num  2 2 2 2 2 2 2 2 2 2 ...
 $ Reach_Num : Factor w/ 27 levels "1","10","11",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Notes     : Factor w/ 74 levels "BoR 10","BoR 12",..: NA NA NA NA NA NA NA NA NA NA ...
 $ geometry  :sfc_POINT of length 5588; first list element:  'XY' num  -113.4 44.7
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA
  ..- attr(*, "names")= chr  "CU_Number" "Fish_Site_" "Site_ID" "StreamName" ...

Structure of my raster layer:

Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots


texts: How to rotate orthographic projection using basemap

texts:
I'm working with the python basemap library and am plotting data over Greenland using an orthographic projection. However I find that due to the location of the coordinates, basemap auto plots the image in this orientation.
Does anyone know how to rotate it so that it is upright instead?
Here is some minimal reproducible code that plots the map without the satellite data:
from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(13, 7))

lon_0 = 42.85565185546875
lat_0 = 75.9834976196289
llx = -2675341.073371827
lly = 316461.22954259627
urx = 286695.76147077326
ury = 3171829.6041646637

m = Basemap(projection='ortho',lon_0=lon_0,lat_0=lat_0,resolution='l',\
        llcrnrx=llx,llcrnry=lly,urcrnrx=urx,urcrnry=ury)

m.drawcoastlines()

plt.show()

Current image:

Desired:


texts: How do I merge two annotations into a single part?


texts: Creating cumulative viewshed (non binary) using ArcGIS Desktop?

texts:
In ArcGIS 10.2.For a viewshed calculation, as an input raster I have a DEM of an area. As an input feature a shapefile with 10 turbines.

Is it possible to create an output raster where the result is NOT just the "visible" or "not visible" one ?

I would like to create a raster that shows how many turbines are visible for each point.


texts: Creating geometry field from lat long fields in Geopackage table

texts:
I am playing around with Geopackages for the first time. I've searched hard for valid syntax to generate a geometry field from lat long fields for an attribute-only table in a Geopackage.
The PostGIS method is
alter table a add column geom geometry(Point,4326);
update a set geom = ST_MakePoint(x,y)
I can't get anything similar to work in the SQLite GPK syntax.


texts: Summarize multiple columns in field calculator

texts:
I am performing the Huff model with 10.000+ facilities and dwellings. As a result, each facility has a value between 0 and 1 that determines the probability a consumer patronages a facility from a certain dwelling. The script generates this value for each store, thus I end up with 10.000+ probability columns. Now I want to add all these probabilities in one accessibility field (per dwelling), since the I am interested in the accessibility from a consumers´ perspective. Is there a function available where all the values are summarized into one field, but exceptions can be given? E.g. Sum all variables horizontally (rows) into one field, where certain columns can be excluded?

Ps. The HUFF-model script also generates a sales potential column per facility in an adjacent column (next to the probability column), which is copied from another cell. This leads towards unnecessary fields. Is there a shortcut to delete all columns that contain facility_X_"SALES". Or can be this be deleted in the HUFF-script?


texts: Openlayers ArcGIS image vs tile service performance

texts:
I want to show my ArcGIS REST MapServer image tiles on my OpenLayers 4 application. But there are two kind of examples to use ArcGIS MapServer.

ImageLayer and ImageArcGISRest
TileLayer and TileArcGISRest

ImageLayer gets large size images like 1390x800 px. But TileLayer gets smaller size images like 256px. I am using dynamic ArcGIS REST service. So I am confused a bit.

Which is the best performance for my client application and ArcGIS Server?


texts: Geopanda write GeoDataFrame into shapefile or spatialite

texts:
I trying to save file from GeoDataFrame to shapefile or spatialite. But for some reason i got error:

My GeoDataFrame appear correctly without any issuesi look's like there is something wrong with my gdf.to_file commend but i have no idea what.
I also checked if i got most up to date lib.

def writeIntoDatabase(input_path, numpy_point_array, 30):

    print ('Reading data into shp')
    df = pd.DataFrame(numpy_point_array)
    dfthresh = df.loc[df[3] < 30]
    geometry = [Point(xyz) for xyz in zip(dfthresh[0], dfthresh[1], dfthresh[2])]
    #df = df.drop(['Lon', 'Lat', 'Alt'], axis=1)

    print(geometry)

    crs = {'init': 'epsg:27700'}

    gdf = GeoDataFrame(dfthresh, crs=crs, geometry = geometry)
    print (gdf)

    gdf.to_file(driver = 'ESRI Shapefile', filename = input_path)

Can you help me?


texts: Android GIS - Reverse Geocode using Shapefile (.shp)

texts:
I want to develop an Android application. I'm still a newbie here.

What I want is to develop my app just like a GIS that can implement reverse geocoding. So I can put in lat / lon coordinates and then text is returned, for example the city name.

But I don't want to implement Reverse Geocode using Geocoder from Google, all I want is for the coordinates that I entered to be checked from a shapefile (.shp) in the database, then return a city name. And I want to use MySQL as a database for my shapefile (shp). Is it possible, and what should I do to make it happen? Could you give a tutorial?

Thanks in advance.


texts: Is it possible to use OSRM in order to calculate routing into a building?

texts: I would to know : If we save in OSRM server a dataset of points stored from a GPS tracker into a building (a flat one) : is it possible to request the server in order to find the shortest way from a point to another ? The aim would be to help customers/employees find their way into a place which has not any standard GIS data like roads, street and so on (like a stadium used during a convention, a zoo, an amusement park, and so on)


texts: Wait for the processing of geometries to be finished in OpenLayer.js

texts:
I need to transform an array of geometries to vector Layers in OpenLayers.
Here is a simple snippet of initial code:
var inputFeatureArray = [geom1, geom2, geom3, geom4, ...]
// I need to populate a second one with new ol.vector.layer
// from the geometries of the first array:
var array2 = [];

var vectorSource = new ol.source.Vector({
    features: new ol.format.GeoJSON().readFeatures(geom1,{
            dataProjection: 'EPSG:4326',
            featureProjection: 'EPSG:3857'
    }),
    projection: 'EPSG:3857'
});

var vectorLayer = new ol.layer.Vector({
    source: vectorSource,
    zIndex: 1,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'rgba(0, 0, 255, 1)',
            width: 2
        })
    }),
});

I was struggling a little before coming to this solution, which sounds good, but is obviously not working.
I have the feeling that the JS code doesn't wait for the second array to be actually filled...
var inputFeatureArray = [geom1, geom2, geom3, geom4, ...]

function processFeatures(inputArray) {
    let outputArray = [];
    for (feat of inputArray) {
        var vectorSource = new ol.source.Vector({
            // This may take some time...:
            features: new ol.format.GeoJSON().readFeatures(feat,{
                    dataProjection: 'EPSG:4326',
                    featureProjection: 'EPSG:3857'
            }),
            projection: 'EPSG:3857'
        });
        var vectorLayer = new ol.layer.Vector({
            source: vectorSource,
            zIndex: 1,
            style: new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: 'rgba(71, 248, 255, 0.9)',
                    width: 2
                }),
                fill: new ol.style.Fill({
                    color: 'rgba(71, 248, 255, 0.2)',
                }),
            }),
        });
        outputArray.push(vectorLayer);
    }
    return outputArray;
}

// Actual call to the function:
var array2 = processFeatures(inputFeatureArray);
console.log("array2: ", array2); // To see what's happening

There is actually no error, but the console.log() is showing empty stuff in my browser:

I don't know how to tell that piece of code to wait for the 2nd array to be actually filled (if that is the problem of course).
In the end, I want to inject all the vector layers into a map:
// Here are the vectors layers in addition to one WMS basemap:
var layers = [tileLayer, array2];

var map = new ol.Map({
    controls: ol.control.defaults().extend([
        new ol.control.ScaleLine({
            units: 'degrees',
        })
    ]),
    layers: layers, // <-- layers plugged here
    target: 'map',
    view: new ol.View({
        projection: 'EPSG:3857',
        center: ol.proj.transform(
            location.coordinates.slice(0,2),
            'EPSG:4326',
            'EPSG:3857'
        ),
        zoom: 8,
    }),
});


texts: QGIS Advanced String Filter Expression returns too many results

texts:
"GNR" is a string column and I expect the filter expression "GNR"='.11' to only return the feature with the value '0.11'.

Why are the other features with values of '.110'and '.1100' included?

Is there something wrong with my query?


texts: Area of polygon calculated as 0 on ArcMap

texts:
I am trying to calculate the area of several polygons on ArcMap (10.3.1). The polygons are all on the same layer.

I am using NAD 1983 UTM Zone 16 as my projected coordinate system. I have tried to calculate the area (m^2) using "calculate geometry..." (my field is set to double with precision 15 and scale 2).
All the areas come back as 0. When I click on individual boxes, I see that the values are not actually 0, but something like 3E-8 meters^2 (which is definitely not correct).

Any suggestions on what may be causing this?


texts: Appropriate map projection for the Pacific Ocean

texts:
I've found this map that shows the Pacific Ocean from http://www.transpacificproject.com:

I'd like to produce a similar map, for the same region of the earth and I'd like to use the same projection as I think it illustrates well the distance between the continents.

I've searched for the name of this projection, but I don't think I know the answer yet. At first I thought it was Robison but according to Wikipedia:

  The Robinson projection is a map projection of a world map which shows the entire world at once.

Which doesn't seem to be the case. (Can't see Europe and Africa here.)

So, my question is two-fold: what is the name of this projection? Should I use it to illustrate what I said earlier or should I use something else?


texts: Auto-updating symbology between *.mxds for same shapefile using arcgisscripting?

texts:
I have a variety of .mxds for any one project I'm working on, which are all different but all use the same set of base shapefiles for my project.

I know I can save a .lyr version of my initial shapefile and import this into subsequent .mxds to retain the original symbology I have set, but if I then change the symbology of the original shapefile, even if I save over the .lyr file with the new symbology, nothing updates in any of my existing files?

My question is, is there an easy way to auto-update this symbology, or is it a manual task? If it's a case of manual re-linking then that's a lot of updating I'm going to have to do on a routine basis - so I'm keen to find a simple way of doing this? I have no python experience but if this is the way to do it then I am eager to learn!!

In terms of what I have found so far for a Python solution:

I found this link here on StackExchange, which I think might be a good starting point in terms of Python, as it seems to be a script to auto-update data sources [ How to update shapefile data source in multiple dataframes and mxds using lyr.replaceDataSource and Python dictionary? ]

And also I found this thread on ArcGIS Idea [ http://ideas.arcgis.com/ideaView?id=087300000008FNBAA2 ] albeit dated 2010, which seems to ask the same question as me, and someone has provided a link to some arcpy.mapping sample scripts which again look very similar to what I need, though they seem to do one .mxd at a time which is about the same effort as relinking them manually.

I did download the sample python scripts to have a look at, and I've installed Python (command line), but I have no experience with this and am struggling a bit with knowing where to start to get these to run!

I'm not sure if there's a simpler way of going about this which AUTO links rather than re-running a script with every edit - it seems odd to me that this is not a function built in to ArcGIS?

I am using ArcMap 9.2


texts:
I'm currently using a combination of ArcMap 10.2 and igraph in R to analyse a network with associated data.  I have an igraph network object based on the OS ITN network with associated TOIDs for each road segment and weights based on distance.  I also have a second network dataset (Shapefile format) that contains TOIDs for each segment, start and end node.

What I'm looking to do is match the second dataset to the first to create a graph object that is weighted according to different criteria (activity counts).  So I need to match the road segments, then create the new graph object, then I'll be able to compare shortest paths (distance vs activity count).

My starting point is something like this...

inds <- match(someIDs, ITNAdjOneWayUpper$IDs)
V(ITNNodeGraphNew)[which(!is.na(inds))]$weights

I'm just not sure what I should be entering to create the new graph object.


texts: Perform affine transformation on shapefile with OGR, GDAL, or other command line tools?

texts:
Transformations with WKT or EPSG code are made easy with ogr2ogr, but are there tools for me to define my own affine transformation to apply to the coordinates of a shapefile?

Or, in general, is there a tool to apply a function to shapefile coordinates?

Further along this note, I also need to change coordinates in the shapefile based on a LUT.

Are there simple tools that exist which do this, or will I have to write such a utility myself?


texts: Running reclassify and Int tools on raster?

texts:
I am trying to do a very simple thing - convert a raster to a polygon where positive values are one set of polygons, and negative values are another set of polygons. The raster ranges from a minimum value of -8.00857 to a maximum of 9.4592.
Here's what I have tried:

Convert floating point Raster to Int (Int does not work).
Multiply raster by a large number to remove decimal places then Convert Raster to Int (Int does not work)
Reclassify into two classes: 0 for negative, 1 for positive (Reclassify does not work)

Extract positive and negative values using Extract by Attribute and then apply steps 1-3 to them. Extract by attribute works, but none of steps 1-3 work.

The error messages do not indicate why the tools are not working.


texts: BBOX and Projections of WMS services in Geoserver

texts:
we have tried to overlay a raster layer on Google map and have some difficulties to align the raster images onto the base map. When it is on the overview (Zoom level 1), it looks like this:

On a lower resolution, you could still spot the mis-alignment in high altitude area such as North Europe and Canada:

BUT when we zoom into the highest level, the problem has gone:

We have tried these approaches to address this problem:

Change the projection of the raster layer on Geoserver from WGS 84 to Web Mercator, does not even shown the image on Google Map.
We try to change the bbox option in the request but it returns strange results that we cannot even zoom in.

Do you have any idea why this happen? What can we do to address this problem?

Thank you very much!


texts: Python OGR: convert a point geometry into polygon

texts:
I'm writing a routine for creating the convex hull envelop around shapefiles of islands.
I could read the data, extract polygons vertices into a geometry and then compute the convex hull.

Now I need to save it to a shapefile of polygons.

Is there a way to transform the geometry (named convexHull in my code) returned by the convexHull method into a polygon geometry; or am I forced to parse all vertices of convexHull and create a ring of a polygon?

def doConvexHull(infile, outfile):
    inH = ogr.Open(infile, 0)
    if inH is None:
        usage("Could not open file {0}. Exit.".format(infile))
    layer = inH.GetLayer()

    # get all polygons
    thisGeometry = ogr.Geometry(ogr.wkbPoint)
    for index in xrange(layer.GetFeatureCount()):
        feature = layer.GetFeature(index)
        geometry = feature.GetGeometryRef()
        ring = geometry.GetGeometryRef(0)
        points = ring.GetPointCount()
        for p in xrange(points):
            lon, lat, z = ring.GetPoint(p)
            thisGeometry.AddPoint(lon, lat)

    convexHull = thisGeometry.ConvexHull()

    drv = ogr.GetDriverByName( "ESRI Shapefile" )
    ds = drv.CreateDataSource( outfile )
    if ds is None:
        usage("Could not create file {0}".format( outfile) )

    lyrname = "convexHull_${0}".format( layer.GetName() )
    lyr = ds.CreateLayer( lyrname, layer.GetSpatialRef(), ogr.wkbPolygon )
    thisFeature = ogr.Feature( layer.GetLayerDefn() )
    thisFeature.SetGeometry( convexHull )
    lyr.CreateFeature( thisFeature ): error from geometry


texts: Learning resources - Getting started with development for ArcGis Online

texts:
I would like to develop search widget / find task for feature layer at ArcGis online website.

I would love to get some resources for that (reading links, videos) for getting started. so far i've install SDK WPF and visual studio, I would like to get resources for better understanding of the process ahead.


texts: Is it possible to have postgres9.6 and postgis 2.3 on ubuntu 14.04?

texts:
I want to upgrade/install postgres9.6 and postgis2.3 on trusty but I cannot find a way...

Has someone succeed on that?

Ps. this repo

sudo sh -c 'echo "deb http://apt.postgresql.org/pub/repos/apt/ $(lsb_release -cs)-pgdg main" > /etc/apt/sources.list.d/pgdg.list'

leads to broken packages issue


texts: How can I view the Certificate of a code-signed ESRI AddIn (*.esriAddIn) without ESRIRegAddIn.exe?

texts:
For users with ArcGIS 10 installed, it's easy to view the certificate for a signed AddIn since it is presented to you in the Confirm Add-In file installation window.

Is there a way to view this certificate information if you don't have ArcGIS 10 installed?

As an alternative, can the ESRIRegAddIn.exe be copied to a computer without ArcGIS 10 and successfully run in order to view the certificate information?


texts: GeoServer Manager: how to change description of stores?

texts:
Here is my code to upload a new layer from a shapefile (zipped):

boolean exists = reader.existsDatastore(workspaceName, storeName);
  if ( !exists ) {
        File shapeFile = new File( "path/to/zipFile.zip" );
        boolean published = publisher.publishShp(workspaceName,
              storeName, "TheLayerName", shapeFile, "EPSG:4326", "");
  }

All is fine, but I need to set the "Description" of the new Store and maybe the "Create spatial index if missing/outdated" value too.

After run this code, I go to GeoServer and select "Stores". On "Edit Vector Data Source" screen I can see the "Description" field is empty.

After googling a lot I discover the GSCoverageEncoder layerEncoder = new GSCoverageEncoder(); method, but I've tried many ways and have no success ( with configureCoverage ).


texts:
This site has excellent resizing (and drifting) animations.
http://www.tomtom.com/livetraffic

Anybody know how the zoom effect is achieved?  OpenLayers has an ancient sandbox with animated zooming, but it was never integrated.


texts: ogr2ogr elasticsearch export to a shapefile

texts:
http://www.gdal.org/drv_elasticsearch.html

Is it possible to export a elasticsearch index into a shapefile using something like:

ogr2ogr -f "ESRI Shapefile" -t_srs crs:84 test.shp ES:http://localhost:9200


texts: I've seen some past comments that a /s option is available to silently install the app. However, that doesn't work with the latest stable version. Does QGIS support a built-in silent install method?


texts: Is there a way to configure QGIS to prevent the st_extent call when adding a PostGIS table?

texts:
I'm trying to add an extremely large table with a points geometry field to QGIS and QGIS appears to hang.  After some investigating in the postgresql database, I tracked it down to a st_extent on everything in the table.

Is it correct to assume that the reason QGIS does this is in order to obtain a full extent for the layer? And if so, is there a way to bypass this? Maybe either use the spatial reference extent or manually provide the extent?

QGIS 2.6 64-bit
POSTGIS="2.1.3 r12547" GEOS="3.4.2-CAPI-1.8.2 r3924" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.0, released 2013/04/24" LIBXML="2.7.8" LIBJSON="UNKNOWN" TOPOLOGY RASTER
PostgreSQL 9.3.5, compiled by Visual C++ build 1600, 64-bit


texts: Wrong/no formatting Word files inserted in ArcMap?

texts:
I'm using ESRI ArcView 9.3.1, and in a project in ArcMap I need to insert a Word file formatted (a table + text, bolded, centered, colored etc.). In Layout View the inserted object shows as in Microsoft Word. But when I'm exporting the map (such as Export to ...jpg, png or pdf) the result with the Word file included is bulk with no formatting.
This error (not retaining the word file formatted) appeared suddenly after no change to the program (ESRI ArcView 9.3.1.) and to the computer (Microsoft Windows 10 64 bits).

I have reinstalled Esri ArcView and the problem isn't solved.

Here is a sample with the problem:

Do you have any idea how to solve this problem?


texts:
I am currently running processing on two feature classes (points, number of points are 5000+ and 8000+ respectively) in ArcGIS 10.2 on personal geodatabase using iteration and selection by attributes/location.

The size of PGDB intially is only about 2MB however after the processing finishes it increases upto 400MB. My outputs of the processing are saved in a separate file GDB.

Can anyone explain to me why the size of original PGDB increases. Is this because the SQL queries are also saved in that PGDB.

How can I delete them? All i need is the original PGDB (with the orginal size) and the outputs stored in a separate GDB.


texts: query.where statement using date range in JS API

texts:
I'm trying to construct a query.where statement using 2 input type=date fields on the page.

The date field in SQL Server in the format 2014-12-31 00:00:00.0000000. The input fields date is in the format 2014-12-31. I am getting errors using the following statement.

function execute () {
    query.geometry = addressPoint;
    startDateQuery = document.getElementById("startDate").value
    endDateQuery = document.getElementById("endDate").value
    query.where = "(DATE_VAL >= " + startDateQuery + ") AND (DATE_VAL <= " + endDateQuery + ")";
    queryTask.execute(query, showResults);
}

If I use the REST interface however,
DATE_VAL >= '01-01-2014' AND DATE_VAL < '12-02-2014'
works fine. So I am not sure if I need to use a formatting option for the query.where statemnent.

Date input fields:

            Start Date: <input type="date" id="startDate" style="color:#333" />

            End Date: &nbsp;<input type="date" id="endDate" style="color:#333" /><button id="SearchLightning" data-dojo-type="dijit/form/Button">Search</button>


texts: Google Earth Engine Matrix Multiply Error?

texts:
In the below code, Tasseled Cap algorithm computed using array based computation but returns error and I don't know how can I solve this:

code link: https://code.earthengine.google.com/7a54500cb06ec88d6ffc6b623965de7c
Map.centerObject(table);
Map.addLayer(table);

var lc8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
.filterDate('2018-01-01','2020-01-01')
.filterBounds(table)
.filter(ee.Filter.calendarRange(7,9,'month'))
.filter(ee.Filter.lt('CLOUD_COVER',10))
.map(function(img){
  var bands = img.select('B[2-7]').multiply(0.0001).clip(table);
  return bands.copyProperties(img,['system:time_start','system:time_end',
  'system:index'])
}).median();

print(lc8);

Map.addLayer(lc8,{bands:['B5','B4','B3']},'lc8',false);

var coefficients = [
  [0.3029, 0.2786, 0.4733, 0.5599, 0.508, 0.1872],
  [-0.2941, -0.243, -0.5424, 0.7276, 0.0713, -0.1608],
  [0.1511, 0.1973, 0.3283, 0.3407, -0.7117, -0.4559],
  [-0.8239, 0.0849, 0.4396, -0.058, 0.2013, -0.2773],
  [-0.3294, 0.0557, 0.1056, 0.1855, -0.4349, 0.8085],
  [0.1079, -0.9023, 0.4119, 0.0575, -0.0259, 0.0252]];

var bands = ['B2','B3','B4','B5','B6','B7'];

var arrayImage1D = lc8.select(bands).toArray();

print(arrayImage1D)

var arrayImage2D = arrayImage1D.toArray(1);

var componentImages = ee.Image(coefficients)
.matrixMultiply(arrayImage2D)
.arrayProject([0])
.arrayFlatten(
  [['brightness','greenness','wetness','forth','fifth','sixth']]);

  print(componentImages);

  Map.addLayer(componentImages,[],'TC',false);


texts: Multiple metadata within a single mosaic dataset

texts:
I am creating a mosaic dataset consisting of thousands of different DEMs.  These will make up a single State, however I need the Counties to be able to hold their own metadata.  Each raster contains its own metadata, however when mosaicked together the mosaic dataset contains no metadata.  I can import an .xml to the description and to the overviews, but this limits me to a single description for each raster within the dataset.  Does anyone know how to make multiple metadata select-able within a single mosaic dataset.

Note: I am using ArcGis 10.1 and sharing these to a ArcGis 10.2 server through an image service.
Each county mosaic works seamlessly, holding its own metadata.  I just cannot figure out how to get different metadata to display when identifying a particular county.

Thanks in advance!


texts: How to generate 512x512 map tiles (2x resolution for retina displays)

texts:
I'm looking for a method to produce 512x512 pixels map tiles using CartoCSS and Mapnik, for example exactly as CardoDB does for its base maps (here below a 256x256 map tile vs its correspondent 512x512 map tile):

I need that for producing 2x resolution map tiles for retina displays. I'm using Openstreetmap Carto map stylesheets and can generate Mapnik xml files.


texts: Counting layer number of raster stack with value under some conditions and write it into output layer?

texts:
I have generated a 8-day resolution gross primary production (GPP) fro several years using MODIS data.Each year has about 46 GPP layers. Now, I want to calculate the carbon uptake periods (in which GPPis larger than 1 g C m-2) based on these data. The result is a single raster layer with the pixel values indicate the number of the raster layers that is larger than 1. I have found a similar answer, but I still cannot solve it.Any response would be help.

The following code is one I used.

GPP2000h26v04=brick(GPP2000049,GPP2000057,GPP2000065,GPP2000073,GPP2000081,GPP2000089,GPP2000097,GPP2000105,GPP2000113,
GPP2000121,GPP2000129,GPP2000137,GPP2000145,GPP2000153,GPP2000161,GPP2000169,GPP2000177,GPP2000185,GPP2000193,GPP2000201,
GPP2000209,GPP2000217,GPP2000225,GPP2000233,GPP2000241,GPP2000249,GPP2000257,GPP2000265,GPP2000273,GPP2000281,GPP2000289,
GPP2000297,GPP2000305,GPP2000313,GPP2000321,GPP2000329,GPP2000337,GPP2000345,GPP2000353,GPP2000361)

CUP.days=calc(GPP2000h26v04,fun=function(x,na.rm) sum(x[x>1]))


texts: Using dynamic text elements from ArcMap in conjunction with text elements derived from Python

texts:
I have a script that cycles through features based on attributes and makes a PDF centered around each feature.  I would like to be able to include both a dynamic text element for the date as well as a title that is derived from the python script.  I'm currently using the  arcpy.mapping.ListLayoutElements(thisMap, "TEXT_ELEMENT") command and then updating the text element for each iteration of the loop.  The problem is that not only does the title get changed but the dynamic date also gets changed to the text meant for the title.  How would I construct the code to only target a specific text element instead of all of the text elements in the map?  I have tried using the wildcard, but the problem with this is that the text gets changed from every previous iteration, so the wildcard will only identify the text from the first iteration.

with arcpy.da.SearchCursor(fc, ("SHAPE@", "UNIQUE_ID")) as searchCur:
for row in searchCur:
        myDF.extent = row[0].extent
        myDF.scale = mapScale

        for elm in arcpy.mapping.ListLayoutElements(thisMap, "TEXT_ELEMENT"):
            elm.text = str(row[1])

        arcpy.RefreshActiveView()
        arcpy.mapping.ExportToPDF(thisMap, str(pdfLoc) + "\Map_" + str(row[1]) + ".pdf", resolution=300, image_quality="Best")   #***test output location variable structure
        arcpy.AddMessage("Output --->  Map_" + str(row[1]) + ".pdf was created successfully!")
        arcpy.AddMessage(str(cycleCount) + "/" + str(recordCount))
        cycleCount += 1


texts: How to export a specific band from an image collection in Google Earth Engine?

texts:
I am trying to export a single band from an image collection as a geotiff file. I am trying following script and each time it is failing.
  // Add the shapefile
  var Tawi = ee.FeatureCollection("users/dsbishtae/Tawi/Boundary");
  var visParams = {color: 'yellow'};
  var nameLayer = 'Tawi';
  var shown = true;
  var opacity = 0.4;

  Map.addLayer(Tawi, visParams, nameLayer, shown, opacity);
  Map.centerObject(Tawi);

  // Load an image **collection**.
  var LC = ee.ImageCollection("COPERNICUS/Landcover/100m/Proba-V/Global")
                              .select('discrete_classification')
                              .filterBounds(Tawi)
                              .map(function(image){return image.clip(Tawi)});

  var visualization = {bands: ['discrete_classification'],
                       min: 0,
                       max: 200,
                       palette: ["282828", "FFBB22", "FFFF4C", "F096FF", "FA0000", "B4B4B4",
                                "F0F0F0", "0032C8",   "0096A0", "FAE6A0", "4E751F", "4E751F",
                                "4E751F", "4E751F", "4E751F", "4E751F","A0DC00", "A0DC00",
                                 "A0DC00", "A0DC00", "A0DC00", "A0DC00","000080"]};

 print(LC)
 Map.addLayer(LC, visualization, 'Landcover');

 Export.image.toDrive({image: LC.select('discrete_classification'),
                       description:'imageToDriveExample',
                       folder:'/gee_files',
                       region:Tawi.geometry().bounds(),
                       scale:100,
                       maxPixels:1e13})

discrete_classification is the band which I want to export. I am aware that the way I am using  Export.image.toDrive is incorrect. I could not get any way to write it correctly. How to do it?
Export.image.toDrive({image: LC.select('discrete_classification'),
                       description:'imageToDriveExample',
                       folder:'/gee_files',
                       region:Tawi.geometry().bounds(),
                       scale:100,
                       maxPixels:1e13})


texts:
I just upgraded to QGIS 2.18.1 Las Palmas.

This comes with a Santa hat on the icon in the top right corner and windows bar.

Is it possible to remove this from my installation?


texts: QGIS with MSSQLSpatial : no geometry. OgrInfo Corrupt data

texts:
I am having this issue while I am trying to connect to my MSSQLSpatial database.
So, I've created this features with just one record with geometry. I can see this feature on Microsoft Management Studio as well on GIS Desktop (Geomedia). My goal is to perform a connection from this database to QGIS. In fact, I can connect all my tables into QGIS  (QGIS version 3.6.0-Noosa) and see their attributes but unfortunately the geometry doesn't show up.
I check the orginfo and you may see the logs as the following:
ogrinfo -al "MSSQL:server=MYSPATIALSERVENAME;database=MYSPATIALDATABASE;tables=MYTABLE;trusted_connection=yes" -fid 1

Result:
Layer name: MYTABLE
Geometry: Unknown (any)
Feature Count: 1
ERROR 1: GetNextRawFeature(): Corrupt data
Layer SRS WKT:
PROJCS["NZGD2000 / New Zealand Transverse Mercator 2000",
    GEOGCS["NZGD2000",
        DATUM["New_Zealand_Geodetic_Datum_2000",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6167"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4167"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",173],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",1600000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","2193"]]
FID Column = ID
Geometry Column = Geometry_SPA
Name: String (50.0)
Geometry: Binary (0.0)
ERROR 1: GetNextRawFeature(): Corrupt data
OGRFeature(testFeature):1
  Name (String) = test1
  Geometry (Binary) = C8FFD20FBC8CCF11ABDE08003601B7696898B09D44DA35418295C131AA06534100000000...

It shows this error: ERROR 1: GetNextRawFeature(): Corrupt data
I also tried with a different table and the same happens. What am I missing here? Why is giving me this corrupt data?


texts: Is it possible to open an attribute table from the Python console or through a Python script in QGIS?


texts: Cannot read pixel value from DTM raster with QGIS

texts:
I am using this DTM:
http://ags.cuzk.cz/arcgis2/services/dmr4g/ImageServer/WMSServer?request=GetCapabilities&service=WMS

At ArcGIS I am able to ready pixel value of the raster but at QGIS it doesn't work. Tried this plugins to get the value

Point sampling tool
Value Tool
Profile tool

But the layer looks empty, like no pixel value in there (at ArcGIS there are).

Can anybody try me link to confirm same issue with layer? Maybe I do something wrong.


texts: Python clipping large shapefile by radius

texts:
I cannot find any real answer so I'd like to ask how to crop shapefile with roads in Python.

The goal is to leave roads (nodes) in black circle only, then output that to dataframe in geopandas (possibly save later as shapefile).

Assuming, that I have some place's coordinates in lat + lon and a very big file of roads, how to clip this file by radius, let's say, of 1km?

Theoretically I can create shapefile of the circle basing on coordinates and then crop roads by this shapefile using some external library like ogr2ogr, but it seems to be very slow and I don't want to use subprocess.

Is it possible to do that directly in geopandas? (I can load all roads to pd dataframe first)

Is there any already implemented option to crop roads but leave all these nodes that are even partially in black circle?


texts: Sampling one shapefile with polygons of another using R?

texts:
I have two shapefiles that partially overlap. One (the bigger one) consists of several polygons which represent fields I am interested in. The second, is a rectangular shape that overlaps about 30% of my first shapefile and hold information about soil texture. I want to clip my first shapefile to the extends of the second and then find out the soil texture for each of the polygons of the first shape.
I tried the merge function and the intersect function in Qgis so far. The merge seemed to be a shapefile consisting only of information of the soil texture but I could not find my fields in it. The intersection gave me an empty shapefile.
Can someone point me in the right direction by providing some procedure names (preferably in R) two achieve what I need. I think I am searching with the wrong keywords here but as a newbie that is what happens;)


texts: MS SQL Layer - Losing LineStrings (Polyline) from SQL Server to QGIS

texts:
I have a table called [gfts].[Shapes] which contains a primary key shape_id and a geography data type GEOM. When I query the table in SQL Server I get the following Spatial Results:

However, when I add an MS SQL Layer in QGIS I only have a few of the LineStrings (PolyLines). Whats the issue here. The data seems to be leaking!


texts: What technology is used behind ArcPy cursors?

texts:
I'm using ArcGIS Server 10.0.
I have a geoprocessing using arcpy search, update and insert cursors.
I'd like to know more about the technology behind these cursors.

Does it use any driver or pool to connect to the database (I'm using Oracle 11g)?

Does it open a connexion to the DB once at the creation of the cursor and close it when the cursor is deleted? What about transactions?

Does the cursor open a new connexion and transaction to the DB each time a row is updated or inserted?


texts: QGIS 3.2 SAGA Raster Values to Points: ERROR: Missing Parameter Value

texts:
When I try to run the saga:rastervaluestopoints tool in the QGIS 3.2 Python console  it returns the following error (see below). I can't figure out what's going wrong. The tool works when running it straight from the processing toolbox.

CODE:

processing.runAndLoadResults("saga:rastervaluestopoints",
                        {"GRIDS": "input.tif",
                        "NODATA ": True,
                        "TYPE": 1,
                        "SHAPES": "output.shp"})

ERROR MESSAGE:

    Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 68, in <module>
  File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\tools\general.py", line 121, in runAndLoadResults
    return Processing.runAlgorithm(alg, parameters=parameters, onFinish=handleAlgorithmResults, feedback=feedback, context=context)
  File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\core\Processing.py", line 133, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Error: Missing parameter value for parameter NODATA        .


texts:
I am currently working on mapping my local fire departments incident responses over the last three years. I have the data in an Excel workbook and would like to bring it in to ArcMap 10.2.2 to analyze and export data as a shapefile. I can get the table loaded to a blank map but when I try to display XY data (right-clicking in the TOC) the system freezes or crashes. When I search the tool box for the make XY event tool it finds it, but when I click on it I receive a script error message:

  An error has occurred in this script on this page. Exception of type
  'System.OutOfMemoryException' was thrown.

What would be causing this error? I am not great with computers but I was wondering if it is because my Excel workbook contains to many features (almost 300,000). Or, is it because my computer is not strong enough to run this data and tools. I have an AMD Athlon(tm) 64 X2 Core Processor 5000+  2.60GHz  3GB RAM  32-bit OS  and running Windows Vista SP2.


texts:
I want to find the least challenging method of getting a shapefile into Excel.

I understand this is often achieved using an .svg file (which must be converted from an .shp).


texts: How zooming Google Maps works regarding the new projection?

texts:
Please help understanding how tiles now work on Google Maps.

As I understand (these are all assumptions, please feel free to correct) the following happens:

Google Maps now work with orthographic projection;
When we zoom from the globe to a territory, the orthographic projection remains; so the first question is, as an orthographic projection is a view from an infinitely distant point of perspective, Google Maps probably had to precalculate a lot of zooms for a lot of territories?
If that's true, what happens when we drug the map? Are all the tiles recalculated? (That's the second question).
The third and the main question is: if the Web Mercator projection is not used anymore, how do tiles zoom? In my understanding, the whole Web Mercator thing was making the map flat, thus cuttable in squares, each square (tile) duplicating each quadrupling each time we zoom. But if the tiles are no longer square, how do they quadruple?


texts:
Im working with MAXENT and im having trouble using some shapefile files with critical variables. I rasterized the shapefiles into rasters selecting one attribute from the attribute table, but before rasterizing, the software asks me to select CELL SIZE. I understand this is the resolution. The world clim .bil that I want to compare to this shapefile has a resolution of: 0,00833 The default CELL SIZE when I try to rasterize the shapefile is 0.15
So:
1) How do I know the shapefile resolution or CELL SIZE?
2) If I change the 0.15 to 0,008333, the raster is generated anyway. How do I know the data I'm obtaining is usable for analysis anyway?

And the other question:
When I rasterize the shapefile, I get a folder with various .adf files

How to I convert these files into a format that the MAXENT can read, such as .grd or .bil ?


texts: Exporting a polygon shapefile to GeoJSON with the longitude and latitude switched to latitude and longitude

texts:
I want to integrate the GeoJSON of a shapefile into React.js to display it as an interactive map with OpenStreetMap used in the background as a tiled image.

I have used QGIS to export the shapefile polygons to GeoJSON format. When opening the GeoJSON file in a text editor the coordinates are arranged as longitude and latitude by default (see below example):

{
"type": "FeatureCollection",
"name": "AB_18_epsg4326",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [{ "type": "Feature", "properties": { "id": null, "BlockID": "18", "AreaName": "Alto Branco" }, "geometry": { "type": "MultiPolygon",
        "coordinates": [ [ [ [ -35.883680840759709, -7.207758761892716 ],
                     [ -35.882709722028402, -7.208009651845575 ],
                     [ -35.882898093859403, -7.208733208791063 ],
                     [ -35.883854465392915, -7.208480923455232 ],
                     [ -35.883680840759709, -7.207758761892716 ]
                ] ] ] } }
]
}

When integrating this to React.js and running the program - the polygon is visualised at the wrong place because the coordinates are supposed to be arranged as Latitude & Longitude in OSM instead of Longitude and Latitude.

Is there way to export this shapefile into GeoJSON with the default positions of the coordinates (i.e. longitude and latitude) changed to Latitude and Longitude using QGIS or R?

The desired export should be:

{
"type": "FeatureCollection",
"name": "AB_18_epsg4326",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [{ "type": "Feature", "properties": { "id": null, "BlockID": "18", "AreaName": "Alto Branco" }, "geometry": { "type": "MultiPolygon",
        "coordinates": [ [ [ [  -7.207758761892716, -35.883680840759709 ],
                     [  -7.208009651845575, -35.882709722028402 ],
                     [  -7.208733208791063, -35.882898093859403 ],
                     [  -7.208480923455232, -35.883854465392915 ],
                     [  -7.207758761892716, -35.883680840759709 ]
                ] ] ] } }
]
}


texts: Get look angle at a specific coordinate on a SAR image

texts:
Can we get the look angle of the SAR sensor on a large image for a given coordinate, if we only have the look angle in the mid-swath (or the near-range & far-range angles)?

Do we get the look angle at a specific location by interpolation or is there a more precise method to estimate it?

If we take as an example of a SAR sensor Sentinel-1, the look angle changes inside the same image from a region to another, especially if the image is large (as is the case of the images taken by Sentinel-1 in the Interferometric Wide Swath mode).


texts: Performing spatial clustering of polygons based on attribute

texts:
I want to perform spatial clustering on the polygon which contains a number of total points inside it. I would like to ask if there is any ways that the polygons can be clustered based on the total number of points inside it. Since I am new to ArcGIS and python, I seek a solution without any code.

My map is as shown below.

The expected results look like this (each color represents a cluster).

But in my case, the cluster should be formed from the regular grids shown.


texts:
I'm using GeoDjango to render many points on a map based on OpenLayers 5.
My points have EPSG 4326 and the map has EPSG 3857, I need to reproject the points.
On the OpenLayers API I've founded two parameters that are useful: dataProjection and featureProjection.

I've used they:

var geojsonObject = {
  "type": "FeatureCollection",
  "crs": {
    "type": "name",
    "properties": {
      "name": "EPSG:4326"
    }
  },
  "features": [{% for d in geometry %}
    {
      "type": "Feature",
      "properties": {
        "pk": "{{ d.pk }}"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [{{ d.coord_lon_lat }}]
      }
    {% if forloop.last %}} {% else %}}, {% endif %}{% endfor %}
  ]
}

var vectorLayer = new ol.layer.Vector({
 source: new ol.source.Vector({
   format: (new ol.format.GeoJSON({
    dataProjection: 'EPSG:4326',
    featureProjection: 'EPSG:3857'
    })).readFeatures(geojsonObject)
 }),
 style: new ol.style.Style({
 ...
 }),
});

var OpenStreetMap = new ol.layer.Tile({
....
});

var map = new ol.Map({
  layers: [
    OpenStreetMap,
    vectorLayer,
  ],
  target: 'map',
    render: 'canvas',
  view: new ol.View({
    center: [0.0, 0.0],
    zoom: 2,
  })
});

But I don't see any points on the map. Why?

If I use the serialization of GeoJSON I can see all points on the map.

var vectorLayer = new ol.layer.Vector({
  source: new ol.source.Vector({
    url: '{% url 'points_geojson' %}',
    format: new ol.format.GeoJSON(),
  }),
 style :new ol.style.Style({
     image: new ol.style.Circle({
         radius: 7,
         stroke: new ol.style.Stroke({
           color: 'rgba(255,0,0,1.0)',
           width: 3,
         }),
         fill: new ol.style.Fill({
           color: 'rgba(200,200,200,1.0)'
         })
     }),
 }),
});

GeoJSON

{
  "type": "FeatureCollection",
  "crs": {
    "type": "name",
    "properties": {
      "name": "EPSG:4326"
    }
  },
  "features": [
    {
      "type": "Feature",
      "properties": {
        "pk": "40"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          14.243133533746002,
          40.840827040200026
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "pk": "41"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          -73.16601634025572,
          40.84706035607121
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {
        "pk": "42"
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          138.602233864367,
          -34.89944783005726
        ]
      }
    }
  ]
}