texts: How to get the exact Circle that user has drawn using Leaflet.draw.Circle?

texts:
I have a simple Leaflet map, where The user can draw a Circle, using Leaflet.Draw.

When you use the following code to get the geometry,

map.on(L.Draw.Event.CREATED, function (event) {
        var layer = event.layer;
        console.log(layer.toGeoJSON());
        drawnItems.addLayer(layer);
    });

You get only the center Point. This issue on Github says that it is expected behavior.
The issue says that you can also get the Radius, and build the Geometry yourself.

I tried to do that using Turf, using

 TurfBuffer(t.circle_layer.toGeoJSON(), t.circleBounds.radius, {
          units: 'meters',
        });

The problem is that the circle created, is different from the Geometry shown on the Leaflet Map. (note the Road intersection on the South-East Side of the circle)

Circle on Screen where user is drawing:

Circle, which is result of Turf buffer:

I guess this is due to the radius being in Web Mercator in leaflt-draw.

I would like to get the circle coordinates that the user has actually drawn. How can I get those?


texts: Planet API not being accepted for quick search but is on stats

texts:
As is similar on another post about this problem I receive the error,

"message": "Please enter your API key, or email and password.",
"errors": []

when trying the script,

#API key stored as an enviroment variable
PLANET_API_KEY = os.getenv('api key here')

item_type = "REOrthoTile"

#API request object
search_request = {
    "interval": "day",
    "item_types": [item_type],
    "filter": combined_filter
}

#first post request
search_result = \
    requests.post(
        'https://api.planet.com/data/v1/quick-search',
        auth=HTTPBasicAuth(PLANET_API_KEY, ''),
        json=search_request)

however I do not get the error when using the script,

#API key stored as an enviroment variable
PLANET_API_KEY = os.getenv('api key here')

item_type = "REOrthoTile"

#API request object
search_request = {
    "interval": "day",
    "item_types": [item_type],
    "filter": combined_filter
}

#first post request
search_result = \
    requests.post(
        'https://api.planet.com/data/v1/stats',
        auth=HTTPBasicAuth(PLANET_API_KEY, ''),
        json=search_request)

another post (planet quick-search api not recognizing authentication?) had this same error and this (https://gis.stackexchange.com/a/239588/139012) response mentioned that a flag may be misconfigured, I cannot find a way to change this though.

To clarify, the API is accepted with the stats search result but not with the quick-search.

I removed my personal api key for security reasons but can send it for tests if needed


texts: QGIS - Continuous buffer image [Heat Map]

texts:
Is it possible to have the representation like the one the heat-map like the figure below using QGIS?
The way seems to use buffer options, but I do not know how.
Does exist any plugin that helps?
I am using QGIS1.8 on windows 7


texts: Select the most recent date from a table using Python/arcpy

texts: Given a dataset with a date field, how would I select the row with the most recent date using Python or ArcPy?


texts: mapnik offline maps using generate_tiles_multiprocess.py

texts:
I've been working for a couple of days to get offline map tiles working for use in JXMapviewer2 but I am having trouble when using generate_tiles_multiprocess.py to process the tiles.

My intention is to allow JXMapviewer2 to process it's background tiles offline for a java project I'm working on. Online tiles work but not on the specific network that it is going to be implemented on.

My progress so far;

Download and install PostgreSQL with PostGIS database.
Update pg_hba.conf to trust
Download osm2pgsql and load pbf file into database using openstreetmap-carto. My tables are; planet_osm* line, nodes, point, polygon, rels, roads, ways and then pointcloud_formats, spatial_ref_sys.

osm2pgsql -k -c -S openstreetmap-carto.style -d gis -U postgres new-zealand-latest.osm.pbf

Edit carto project.mml to contain PostGIS connection info
Download and place shapefiles in correct carto directory
Generate XML file and remove font references

carto project.mml > osm.xml

Update generate_tiles_multiprocess.py to point to my osm.xml file and tile_dir. Include bbox and zoom levels.

if __name__ == "__main__":

    #home = os.environ['$HOME']
    try:
        mapfile = os.environ['MAPNIK_MAP_FILE']
    except KeyError:
           mapfile = 'C:/Mapnik Tiles/osm_tiles/openstreetmap-carto-master/osm.xml'
    try:
        tile_dir = os.environ['MAPNIK_TILE_DIR']
    except KeyError:
        tile_dir = 'C:/Mapnik Tiles/osm_tiles/openstreetmap-carto-master/tiles/'

    if not tile_dir.endswith('/'):
        tile_dir = tile_dir + '/'

    #-------------------------------------------------------------------------
    # Change the following for different bounding boxes and zoom levels
    #
    # Start with an overview
    # NZ
    bbox = (166.245, -47.339, 178.660, -33.943)
    render_tiles(bbox, mapfile, tile_dir, 0, 18, "World")

Run generate_tiles_multiprocess.py

This is where I have the issues. The folder structure gets created, eg. 0, 1, 2 etc but no tiles are created.
The console output displays multiple lines of "* attribute 'repeat-distance' with value '50' at line 39259" and sits. (more than 999 lines of this with varying line references)
I then have to Keyboard interrupt to get back to command line.

If I use generate_image.py it creates a map of New Zealand and shows the same attribute lines as above.

Is it a case of leaving it to create the tiles?

I've left it for 10 minutes with no images being output.

JXMapviewer2 java example code;

// Create a zip-file based TileFactoryInfo for OpenStreetMap
// See https://docs.oracle.com/javase/8/docs/api/java/net/JarURLConnection.html for URL syntax details
// You can create such a ZIP file by zipping the content of ${HOME}/.jxmapviewer2/tile.openstreetmap.org
// The ZIP file should contain a list of folders (0, 1, 2, 3 ...) that represent the OSM path structure

TileFactoryInfo info = new OSMTileFactoryInfo("ZIP archive", "C:/Mapnik Tiles/osm_tiles/openstreetmap-carto-master/tiles/tiles.zip!");


texts:
I am hoping to come up with a script that would allow me to mosaic imagery into a polygon tile grid.  My hopes is to do this with gdal or QGIS.  I am wondering if anyone knows of an existing script, before I start.  My first thoughts was to use gdal_merge, but I do not see a variable for using a polygon (shapefile) as the boundary.  gdal_merge does have the -ul_lr function, so I know that it will work on subsets of images.

What I am looking to do is the following: From a source of ortho imagery, transfer any and all pixels that intersect with the tile grid, then move on to the next grid.  There may be 10 source images (for example) that are large in size (geographic extent and disks space) that would be used to populate the individual grid tiles (ie: 1km X 1km).

Thanks....


texts: Projecting Data with Make XY Event Layer?

texts:
I have a directory full of .CSVs representing point data. I run this code on them:

#define output feature class path
outFC2016 = os.path.join(r"...\Output.gdb\outputFDS", "Data_2016")
if arcpy.Exists(outFC2016):
    arcpy.Delete_management(outFC2016)

#define table directory
table16dir = r"...\Data_2016_Dir"

#list all tables in dir
tables16 = os.listdir(table16dir)

#define MakeXYEventLayer_management parameters
in_x_field = "Easting"
in_y_field = "Northing"
in_z_field = "Altitude"
sr =  r"Coordinate Systems\Projected Coordinate Systems\Utm\Nad 1983\NAD 1983 UTM Zone 11N.prj"

#create empty list to hold lyrs for merge
inTableLyrs16 = []

for table in tables16:
    if table.endswith(".csv"):
        inTable = os.path.join(table16dir, table)
        lyrName = table.replace(".csv", "") + "_lyr"
        if arcpy.Exists(lyrName):
                    arcpy.Delete_management(lyrName)
        inTableLyrs16.append(lyrName)

        #make xy event
        arcpy.MakeXYEventLayer_management(inTable, in_x_field, in_y_field, lyrName, sr, in_z_field)

#merge all 2016 GOEA data
arcpy.Merge_management (inTableLyrs16, outFC2016)

It produces a feature class that draws near the equator, when it should be drawing in Southern California.

I think that the spatial reference is correct for a few reasons. First, the metadata for the table says the coordinates are in NAD83 zone 11. That should be 11N. Second, the value I have for "sr" is taken directly from the ESRI example online. Third, I'm pushing it into a feature dataset with UTM 11N projection. Fourth, when I use the Go To XY tool in ArcMap, it places any point from one of these tables in the correct location. Fifth, if I use Add XY Data, and assign the projection to UTM Zone 11N, they show up in the correct location.

One possible factor worth mentioning is I'm working on Citrix.

Can anyone tell me why it doesn't work when scripted?

Here is a sample of the table I'm working with.

Easting     Northing    Zone    Altitude
521341.6335 3618762.765 11      901
521396.3151 3618613.56  11      917
522064.6776 3618697.498 11      975
522063.268  3618699.158 11      981
522062.5182 3618699.156 11      980
522066.8389 3618694.953 11      970
522064.0168 3618699.603 11      977
522063.2746 3618696.165 11      981
522064.6766 3618697.942 11      980
522066.0844 3618697.058 11      975
522004.2594 3618554.804 11      1069
522004.2678 3618551.035 11      1068
522005.6709 3618552.257 11      1058
522004.9238 3618551.036 11      1068


texts: How to dissolve polygons that have been split back into their original polygons?

texts:
I had a feature class containging about a thousand polygons.After some processing such as intersecting, dissolving, merging, more intersecting, I ended up with a feature class of about 6000 polys. This was several months ago, and since then more data has been added/removed, attributes have changed etc.
I now need to get the feature class back to its original state. Not original as in the actual shapes from before, but where one original polygon was split into 4 polys, it must now be dissolved back into one.

One method I considered was using Spatial Join with the match option of Shares_a_line_segment_with, as all the polys which share a line segment would have originated from the same polygon. Unfortunately, the unique IDs assigned to the original polygons were 'lost' a while ago, so I cannot just dissolve the result of the spatial join based on target or join ID, as multiple polys share line segments, and they would all be joined to each other (resulting in several join and target IDs in overlapping polys).
I don't need to preserve the attribute information, I only want the original polygons back. I was thinking of using a feature selection iterator to select by location on each feature where it shares a line segment with another, populating a list with the relevant IDs and then calculating a unique ID for the group, which I could use as my dissolve field, but that seems unnecessarily complicated, and I'm sure there is something simpler I could do.


texts: What does Unit-Scale actually do in GEE?

texts:
I am experimenting with Google Earth Engine and normalized my images using UnitScale as per this post.

When I view the image in the lower pane in Google Earth Engine Code editor and switch the visualization parameters to 100% (from custom: 0-1 - set by the code), I notice that the min and max are nolonger 0-1 but some other range such as 0.03 - 0.56.  Sometimes as I pan and zoom, the min and max of 100% changes.

I need to be certain of what is happening here so I can know whether I can trust mathematics on the normalized data. So my question is this:

If your data range already falls between 0-1 (ie min>0 and max<1), does any scaling take place?  It is not clear from the documentation for this use-case.
Why do the min and max values change when panning and zooming?  I assume this is because the min and max in the pop-up visualization parameters only refere to the portion of the image that is on-screen at a given time - which means that what you see is a 'moveable feast'  - zooming out to 100% does not give a range of 0-1 so I assume this is due to resampling at low zoom levels.  Are these assumptions correct?


texts: Generate random points in geometry near heatspots

texts:
I have an algorithm, that generates random points with in a given polygon (Generate random points in polygon. Algorithm question). The Polygons the points are generated in are communities or districts of germany. So instead of just randomly placing the points somewhere in the Polygon I would like to place them near cities, town etc. since a Point should represent a person.

So my Idea was to use OpenStreetMap data to query cities, town etc., but I yet have to find out which data from OSM I need to use. Under the assumption that the data is available how do I influence the algorithm to generate the Random points near the cities?


texts: How to use QgsLayerTreeViewDefaultActions Class

texts:
I am writing plugin, that connects to postgresql and loads some layers.
I am trying to add my own QgsLayerTreeViewMenuProvider for some groups with my layers but to let old menus for all other layers and groups.

I can not understand how to use QgsLayerTreeViewDefaultActions. Can anybody show some example?


texts: Filtering anomalous pixels from masked raster in Python?

texts:
I am struggling to succesfully filter out anomalous pixels from a Sentinel-1 image band mask (Pixel based), produced in Python.

I am new to Python.

By anomalies, in this situation, I mean individual (or small clumps of) pixels that give off the same backscatter value as water features, but are most likely something else.

I'm currently working on a large scale water body mapping project using SAR (Sentinel-1), using the backscatter as the identifier of permanent water bodies.
The pre-processing has all been done using SNAP, using batch processing, but my mask has been created using Python (GDAL and Numpy), due to the amount of images I will be doing this for.

This is my code:

from osgeo import gdal
import numpy as np
ds = gdal.Open(r"Processed_Sentinel1.tif")
ulx, xres, xskew, uly, yskew, yres =ds.GetGeoTransform()
band_ds = ds.GetRasterBand(1)
Sentinel1 = band_ds.ReadAsArray()

Mask_Zero = Sentinel1 > -21  #In units of the raster, e.g. Pixel Value
Ones_array = np.ones((Sentinel1.shape[0],Sentinel1.shape[1]),dtype=np.float32)

Ones_array[Mask_Zero] = 0

out_raster = gdal.GetDriverByName("GTIFF").Create(r"Sentinel1_Mask.tif",Sentinel1.shape[1],Sentinel1.shape[0],1,gdal.GDT_Float32)
out_raster.SetGeoTransform((ulx, xres, xskew, uly, yskew, yres))
band = out_raster.GetRasterBand(1)
band.WriteArray(Ones_array)
out_raster = None
band = None

The following code has produced reasonable results, but with a fair few "false positives" or "anomalies", and for only one image at a time - Please see below:

I know water isn't the only surface that gives off a low backscatter, hence the results, but it would be useful if I could filter out pixels and produce a smoother map. For example, in SNAP, I would apply a median 5x5 filter, and this would usually remove the bulk of the "False Positives".

How would I filter out these "False Positive" pixels according to their surrounding pixels, in Python?


texts: Not Able to Pass Where Clause in SearchCursor in ArcPy

texts:
Can you please take a look at this snippet and let me know why I am not able to properly pass the Where Clause ('"[NAME_1]" = Ohio') in SearchCursor?

import arcpy
from arcpy import env
def unique_values(table , field):
 with arcpy.da.SearchCursor(table, [field], '"[NAME_1]" = Ohio') as cursor:
 return sorted({row[0] for row in cursor})
uniques = unique_values(r'C:\arcgis\ArcTutor\AAA\src\USA.shp' , 'NAME_2')
for unique in uniques:
  print (unique)

I am getting this error on runtime

Traceback (most recent call last):
  File "<module1>", line 8, in <module>
  File "<module1>", line 7, in unique_values
  File "<module1>", line 7, in <setcomp>
RuntimeError: Unspecified error

Same error on this format as well

with arcpy.da.SearchCursor(table, [field], '"NAME_1" = Ohio') as cursor:


texts: How to save and reload legend settings using Python?

texts:
I know how to change settings of a legend individually by arcpy (like: Patch, Title, Items, Columns, Frame, Element Size&Position, etc.), but now I would like to use the exact legend settings another time and load them using ArcPy.

Is there a simple way to save/load these legend settings all together in/from a text file?


texts:
Pretty straightforward question here.

At an increasing rate I see crowd sourcing data being mapped to great effect.  Where are the raw data sources for all these maps?  I know they vary but I don't think I have even found a single source.

For instance, twitter maps, where is the data?


texts: When i click on Geoserver Web Admin Page gives Message as 'Bad Request (Invalid Hostname)' in web page. I have already data of on going project in geoserver. I have started Geoserver. My geoserver version 2.1.3. in window Xp operating System. Can any one please guide me to get out of the error and get back Geoserver web Admin Page. Please see the Page result in images below.


texts: How to create a QGIS processing script that adds a sequence to a unique identifier column in PostGIS?

texts:
Can anybody help me to create a QGIS processing script that adds a sequence to an existing unique identifier column (type: integer) in PostGIS?

This would be quite helpful, e.g. as a workaround for bug #6798. Unfortunately, I do not have any Python experience.

CREATE SEQUENCE /*input_schema*/./*input_table*/_/*uic*/_seq OWNED BY /*input_schema*/./*input_table*/./*uic*/;
SELECT SETVAL('/*input_schema*/./*input_table*/_/*uic*/_seq', (SELECT MAX(/*uic*/) FROM /*input_schema*/./*input_table*/));
ALTER TABLE /*input_schema*/./*input_table*/
ALTER COLUMN /*uic*/ SET DEFAULT nextval('/*input_schema*/./*input_table*/_/*uic*/_seq'::regclass);


texts: Exporting High Quality TIFF Files on ArcMap

texts:
I am currently working in ArcMap 10.5.1 and am wanting to export a Map as a TIFF file with increased resolution.  If I were to increase the quality of the image to anything above 96 DPI then the image will just appear blank after exporting.

Here are my current steps in ArcMap 10.5.1

Here are my current steps:

1)  File (On Upper left tool bar) => Export Map

2)  Chose my folder that I want the image saved in

3)  In the Options section.  I increase the Resolution by increasing the DPI

4)  Then I press save and the image is in my folder

The above steps only work if I leave the resolution at the default DPI (96 DPI).  If I were to increase that number, after the image is saved, it will just appear blank.  Do you know why this is happening? And do you know how I can create a high quality tiff file in ArcMap?


texts:
I have a problem using a Web Map Tiling Services in QGIS.
Here you can see the address of the service: https://wmts.geo.admin.ch/1.0.0/WMTSCapabilities.xml?lang=fr
I can easily add it through the "Data Source Manager" and I have all my different tilesets. If I load one of them from the "Data Source Manager", a window appears where I need to set a dimension (time).

It works perfectly...
But now, if I go in my browser panel, I also have all my tilesets.

But if I load a data from the browser panel, the layer is empty?
If we look the source in the properties of the layers, we can see that they are also different:

The one who works:

The one who doesn't work:

Did anyone have an idea, so I can load my data from my browser panel? It is boring to go in the "Data Source Manager"...


texts: Issue with gdal getprojection and setprojection returning different information

texts:
I am using Landsat imagery from the USGS. Those who use them may know that the projection attribute is wrong (everything is in the Northern hemisphere).
I have a script to process tar.gz files and output a stacklayer with selected bands of interest. I deal with the projection error by modifying the EPSG value and writing a new projection like this:

orgeoproj =  dst.GetProjection() # North projection
osn = osr.SpatialReference()
osn.ImportFromWkt(dst.GetProjection())
nEPSG= int(osn.GetAttrValue("AUTHORITY",1))
sEPSG= nEPSG + 100
print "mudando EPSG Norte a Sur"
osn.SetFromUserInput("EPSG:%d"%sEPSG)
modgeoproj = osn.ExportToWkt()

When I print orgeoproj (original geoproject), the output has Authority values like:

PROJCS["WGS 84 / UTM zone 23N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32623"]]

After modifying the EPSG value, the modgeoproj is OK now and showing the right EPSG value for the southern hemisphere! The problem rises when I write the modgeoproj with SetProjection(modgeoproj) to the stacklayer and then load it with another script.

After a GetProjection(), I only get:

PROJCS["WGS_1984_UTM_Zone_23S",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]

There's no Authority Information and I need it!
What can possibly be wrong?


texts: Using JavaScript and cartodb.js to get layer's total amount of visible data in certain moment?

texts: I am working with a torque layer and I would like to know for each moment the amount of data points shown in the visualization.


texts: Permanently setting polygon/polyline segment length using ArcGIS Desktop?

texts:
I'm trying to create a several hundred polygons in ArcGIS 10, each of which may have up to several hundred vertices. For consistency sake I would like each line segment in these polygons to be of a set length, say 10m.

I know you can right click and set the next segment length (or use the keyboard shortcut Ctrl+L) as you are drawing, but this is going to insanely repetitive if I have to do this thousands of times. Is there a way to permanently set the segment length as you are drawing a polygon or polyline?


texts: How do I change polygon opacity using a esri WFS in leaflet?

texts:
Note...Solved and here for anyone having this problem...damn simple problem, I had not used a higher case character for fillcolor and fillopacity....they should be fillColor and fillOpacity  , it works after changing that.  mmm

I seem to be able to change the opacity of the polygon outline but not the fill itself.  Can find what is going wrong with my code.  I try to use the fillcolor and fillopacity but nothing changes on my map?

var LMDB1 = new L.FeatureGroup();

    L.esri.featureLayer("https://services.arcgis.com/HIDDEN/arcgis/rest/services/HIDDEN/FeatureServer/1", {

        style: function (feature) {
        if(feature.properties.TAG === 'Lowest Health'){
            return { color: "#AD0909", weight: 1, opacity: .1, fillcolor: "#000000", fillopacity: 1};
            } else if(feature.properties.TAG === 'Low to Average Health'){
            return { color: "#FF8000", weight: 1, opacity: .1 };
            } else if(feature.properties.TAG === 'Average Health'){
            return { color: "#FFFB0D", weight: 1, opacity: 1 };
            } else if(feature.properties.TAG === 'Average to High Health'){
            return { color: "#C4FFD3", weight: 1, opacity: 1 };
            } else if(feature.properties.TAG === 'Highest Health'){
            return { color: "#0E8F2E", weight: 1, opacity: 1 };
            } else {
            return { color: "#FFFFFF", weight: 1, opacity: 1 };
            }
        },
        onEachFeature: function (feature, layer) {
            layer.bindPopup(L.Util.template(poly_templ, feature.properties));}

        }).addTo(LMDB1);


texts: Writing 2 for loops by .map() function in GEE

texts:
I would like to calculate the annual Hansen loss for different elevation classes.
I want to write a function that gives me a list holding 18 years of loss per elevation class.
I have written the following lines:
var loss = function(i) {
  var losscalculate = ee.Image(i) //for each year
    .updateMask(elclasses.get(k)) //for each elevation class
    .reduceRegion({
      reducer: ee.Reducer.sum(),
      geometry: peru,
      scale:30,
      maxPixels:1e9,
      bestEffort: true });
    return losscalculate ;
};

Here I am struggling with writing a proper .map() line. Is there any way to use two variables in a .map() function. Or do you have any suggestion for me to write a 2 loops function?


texts: Uploading to ArcGIS Online using the Cache toolset

texts:
I am using the following steps in ArcMap:

Generate Tiling Scheme tool
Manage tiling scheme tool
Export Tile tool

in order to create a tile package for upload to ArcGIS Online.

I know my map document is more than 200MB( a large dataset), however, only 4MB is uploaded to ArcGIS Online.  I always realize the file is 4MB after I use the Manage tool.

Can anyone tell me why this is happening?


texts: Centroid coordinates for odd-shaped polygons

texts:
I have a problem in computing the centroids of certain polygons in QGIS (via the geometry tools): for certain polygons - with a strange and rather long shape - the centroids computed are off the polygon. Look, for instance, at the polygon marked in red, in the picture below:

Its centroid (pointed at by the blue arrow) is actually falling outside the polygon of interest. This occurs for a handful of (yet annoying) other polygons.

I am relatively new to QGIS, so I don't think whether this problem is due to the projection of the shapefile with the polygons (WGS 84). I tried anyway by re-projecting the map with WGS84/UTM Zone 31, which corresponds to Belgium (the country I'm looking at).

I have not found an answer to a similar questions in past discussions; apologies if I accidentally overlooked it.


texts: Extract coordinate values from Leaflet marker to Java Script variable

texts:
I am trying to incorporate pgRouting layer from geoserver to my HTML GUI. For this I have been trying to make a function that will take coordinate values from a Leaflet marker and feed them into my specified parameters. The parameters will act as start and end points for the routing layer. I can't seem to get any lat/lon values to enter into my var viewparams. The code is as follows:
var params = {
    LAYERS: 'IST_Mosaic:pgRouting',
    FORMAT: 'image/png',
    SERVICE: 'WMS',
    request: 'GetMap'
}
var selectedPoint = null;

var sourceMarker = L.marker([33.5183, 73.1789], {
    draggable: true
})
.on("dragend", function(e) {
    selectedPoint = e.target.getLatLng(); console.log(selectedPoint)
    getRoute();
})
.addTo(map);

// draggbale marker for destination point.Note the marker is initialized with an initial destination positon
var targetMarker = L.marker([33.5191, 73.1768], {
    draggable: true
})
.on("dragend", function(e) {
    selectedPoint = e.target.getLatLng(); console.log(selectedPoint)
    getRoute();
})
.addTo(map);

function getRoute(){
    var start = sourceMarker.getLatLng();
    var end = targetMarker.getLatLng();
    var viewparams = [
        'x1:' + start[0], 'y1:' + start[1],
        'x2:' + end[0], 'y2:' + end[1]
      ];
      params.viewparams = viewparams.join(';'); console.log(viewparams)
      var routing = L.tileLayer.wms('http://localhost:8010/geoserver/wms' , {
        layers: 'IST_Mosaic:pgRouting',
        format: 'image/png',
        transparent: true
    });
}
getRoute();

Please can anyone suggest a solution to this and point out what I am doing wrong. A screenshot for the results:


texts: QGIS 2.10.1 Copy and Paste Problem - Attributes not showing

texts:
I have two shp layers with different numbers of fields and slightly different attributes of those field, like the column name or field length.
I arranged them so they are in the same order and now want to copy points from the Layer Georef_FO to the Layer Georef_FO_2.
Screenshot attached.

I put both layers in edit mode while in the attribute table window and then just copied all records from Georef_FO and pasted them into Georef_FO_2.
I get all the points on the map but no information in the Name Field, just NULL Values.
(the "Name" colums is right now the only column in Georef_FO containing any information and I want this to be copied into the layer "Georef_FO_2".

Restarting QGIS doesn´t make the values appear in Georef_FO_2. There are still only NULL Values.

I couldnt find out how to change field length of a field in an existing layer, so if that may be helpful please explain how to do that.


texts: Problem – QGIS 2.2 - Values under Color Rendering saved in style not applied when applying the loaded style?

texts:
I have saved two variants of the same color-scheme by putting in different values under “color rendering” in the lower part of the screen. It’s a singleband pseudocolor scheme.

Now when I load each of the styles as qml-File while in the Properties-Menu of my DEM it instantly shows me the desired colors in the canvas. I can also see the slight difference between the two variants in the canvas resulting from the different “brightness” values in my two styles.

In this "Preview Mode" everything is great.
But the Moment I click on "Apply" or "OK" when I want to finally apply one of the two versions it switches back to a state as if no changes under color rendering have been made.

So is this a bug? When first loading the style it shows the correct adjustments (so it must be stored somewhere inside the qml file) but when clicking "Apply" all Information under color rendering is ignored?

Or am I doing sth wrong?

When I open the properties windows again on the quite dully colored DEM after trying to apply my style and then manually reenter my desired value under "color rendering" and click "Apply" it thens also applies the desired tweaks.

Edit:
With a Singleband Grey Style it works. I have two variants with different values under "color rendering" and when loading and applying everything looks correctly.


texts: Raster grid on PostgreSQL 12.2/PostGIS 3.0.1 not able to add to QGIS 3.12

texts:
I'm trying to do the example on PostGIS Cookbook by uploading the Temperature Grid exactly as shown in the book (https://subscription.packtpub.com/book/big_data_and_business_intelligence/9781849518666/1/ch01lvl1sec16/importing-raster-data-with-the-raster2pgsql-postgis-command)

$ raster2pgsql -I -C -F -t 100x100 -s 4326 worldclim/wc2.1_10m_tmax_01.tif tmax01 > tmax01.sql
$ psql -d postgis_cookbook -U me -f tmax01.sql

The upload seems to be successful (I'm using PostgreSQL 12.2/PostGIS 3.0.1) on Windows 10 platform (Localhost).

However, when trying to add the resulting image from the database into QGIS 3.12, I get the following error on top of the viewer window:

Invalid Layer: Raster layer Cannot instantiate the 'postgresraster' data provider

Any clues?

Regards,


texts: What does n-D stand for in the following operator definitions?

texts:
PostGIS shows a listing of operators that make reference to n-D

&&& — Returns TRUE if A's n-D bounding box intersects B's n-D bounding box.
&&&(gidx,geometry) — Returns TRUE if a n-D float precision bounding box (GIDX) intersects a geometry's (cached) n-D bounding box.
&&&(gidx,gidx) — Returns TRUE if two n-D float precision bounding boxes (GIDX) intersect each other.
<<->> — Returns the n-D distance between the centroids of A and B bounding boxes.
<<#>> — Returns the n-D distance between A and B bounding boxes.

What does "n-D" mean?


texts: How to select a feature with the mouse click programmaticaly

texts:
I am working in c++/QT 4.8.5 and the api c++ qgis 2.8.2

In my application i had top select layer with the mouse click button to get attribute value and show them in a popup windows

I see in the api online documentation taht there is two classes :

QgsMapToolIdentifyFeature and  QgsMapToolIdentify

can i use one of them to do what i want to do in my application or there is another to do it ?


texts: Create a pivot table for m-n relationship in QGIS?

texts:
My dataset comprises many polygons representing workspaces over a large area:

-1 vector layer comprising my polygons (each polygon ID appears 1 time)

For my online atlas however I need to obtain a pivot table layer as follows :

-1 pivot table listing all the polygon couples less than 300 meters apart (each polygon ID appears over more than 1 line depending the nb neighbor polygons).
I have no clue how to obtain this in QGIS.


texts: Style features based on difference between current date and date in attribute table

texts:
In QGIS >2.18, is it possible to style features based on the age gap between the current date and date recorded in the attribute table. For instance I have purchased a significant amount of OS data for different mapping projects with an annual licence. I have a shapefile which hosts the extent of the purchase area and includes the expiry date for each licence/feature. I would like to apply a graduated style to the shapefile which will give me a visual indication if the licence for a particular area has expired or not. The shapefile can be loaded into QGIS from a spatialite database if it is necessary.

Just wondered if QGIS has the capability to run an expression based off the current date and update the styling every time the layer is loaded into the canvas.


texts:
Is there a way to get line and polygon data without the MultiGeometry Wrapper from a geoserver layer when accessing via the KML Reflector method? It always generates additional Point data for Lines and Polygons. QGIS unfortunately doesn't support wrapped MultiGeometry data. So in QGIS this Data is not valid (Error Message: Unsupported geometry type 'Unknown type' in layer ...):

<MultiGeometry>
      <Point>
         <coordinates>15.007516282814333,48.86425868847214</coordinates>
      </Point>
      <LineString>
         <coordinates>15.00746377393497,48.86434478906123 15.00748822834979,48.86430296174839 15.00757027103822,48.86417348864424</coordinates>
      </LineString>
</MultiGeometry>

While this would work:

<LineString>
     <coordinates>15.00746377393497,48.86434478906123 15.00748822834979,48.86430296174839 15.00757027103822,48.86417348864424</coordinates>
</LineString>

The data is stored in a postgres DB. The geometry type in the sql query for the geoserver layer is set to LineString and EPSG:4326. The sld used on this layer is the default simple blue line.

Any ideas how to access KML Data from geoserver so it can be read by QGIS? Running and converting the data using ogr2ogr is not an option (the data must be accessible via REST).


texts:
I am new in sql and i am trying to get that works a subquery with pgr_drivingdistance and pgr_pointsaspolygon but i get the next error:

  subqueries in FROM must to have an alias,for example, FROM (SELECT ...) [AS] foo.But in the code below the alias is di

Here is the code:

SELECT 1 as result,* FROM  pgr_pointsaspolygon('SELECT id, ST_X(the_geom) AS x, ST_Y(the_geom) AS y FROM (
     SELECT id,seq, id1 AS node, id2 AS edge, cost, the_geom FROM pgr_drivingdistance(
       ''''SELECT id, source, target, time_4 as cost FROM jaen'''',
          1, 15, false, false
      )as di
      JOIN jaen_vertices_pgr pt
      ON di.id1 = pt.id)');


texts: Combining 2 raster data sets with different extent, resolution and point regularity

texts:
I'm having trouble merging two raster datasets that have different resolution and extent, and have irregular point data. Below is info on each raster. I need these datasets to merge according to the grid points so that I can run a fire spread model (which requires both wind and slope data).
I have tried converting to normal dataframes (using rasterToPoints) before merging but the differences lead to loss of a lot of grid points.
I have tried to align the rasters with project raster and rasterize, but I haven't managed to get it to work.
Data from https://globalwindatlas.info/downloads/gis-files
> wind <- raster("ZAF_wind-speed_10m.tif"); wind
class      : RasterLayer
dimensions : 11271, 11804, 133042884  (nrow, ncol, ncell)
resolution : 0.0025, 0.0025  (x, y)
extent     : 13.33407, 42.84407, -50.31423, -22.13673  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : C:/Users/s2000128/Documents/Flammability modelling/ZAF/ZAF_wind-speed_10m.tif
names      : ZAF_wind.speed_10m

Data from https://datacatalog.worldbank.org/dataset/world-slope-model
> slope <- raster("ZAF1_msk_alt.grd"); slope
class      : RasterLayer
dimensions : 1548, 1992, 3083616  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 16.4, 33, -34.9, -22  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=WGS84 +no_defs
source     : C:/Users/s2000128/Documents/Flammability modelling/slope_deg_0/slope_deg/ZAF1_msk_alt.grd
names      : ZAF1_msk_alt
values     : -26, 3264  (min, max)

Here are some examples of what I have tried:
1
windresampled <- projectRaster(wind,slope,method = 'ngb'); windresampled

2
wind_points <- rasterToPoints(wind)
coordinates(wind_points) = ~x+y
proj4string(wind_points) = CRS("+init=epsg:4326")
gridded(wind_points) = TRUE
g_wind <- raster(wind_points); g_wind
extent(wind) <- c(16.4,33, -34.9,-22)
res(wind) <- c(0.0025, 0.0025)
r_wind <- rasterize(wind_points, g_wind, field = wind_points$z, fun = mean, na.rm = TRUE)


texts: Raster calculators in SAGA return incorrect results

texts:
Here's a sample of my dataset (one cell) with 28 different raster values, representing 28 individual Float32 rasters with identical parameters.

0.0000264434,   0.0049953502,   0.0007529828,   0.0002075955,   0.0000516699972,    0.0000000000,   0.0053563155,   0.0000000000,   0.0000000000,   0.0073998945,   0.0000804775,   0.0000000000,   0.0000387653,   0.0000000000,   0.0000000000,   0.0000000000,   0.0000000000,   0.0000000000,   0.0007369027,   0.0000519922,   0.0053393530,   0.0000000000000,    0.0000000000000,    0.0000000000000,    0.0000808992875,    0.0000000000000,    0.0000000000000,    0.0000000000000,

Due to system overload when I run the algorithm in a modeler, I've split the process into 4 groups of 7 layers by order of appearence.
Each subgroup does the summary correct in all the calculation algorithms, including both SAGA ones: Raster Calculator and Grids Sum.
However, when I do the final summarization of the 4 rasters, It appears that the SAGA algorithms return a different result

Excel: 0.0251186417
GDAL_CALC: 0.0251186417
SAGA Raster Calculator: 0.060915395617485046
SAGA Grids Sum:0.0012134892167523503

It's worth noting that if if I run both SAGA algorithms in a standalone way on all 28 rasters, I get the correct result like excel, gdal.

I tried to look if the SAGA algorithms truncate the data somehow but it doesn't seem to be the problem.

Any ideas why this is happening?


texts: Change OpenLayers 3 script into OpenLayers 2.13

texts:
I want to change the code from pgrouting latest workshop which is used the OpenLayers 3 beta to OpenLayers 2.13. The code is in here http://workshop.pgrouting.org/chapters/ol3_client.html

I tried to change it, but still have got some errors like this one: startpoint.getPosition is not a function

Is there anyone could help me? Because there are some features in OpenLayers 2.13 which are not (yet) supported in the Ol3js. My idea is simple, I just want to have a Layerswitcher button, and adds radio button for start and final location for the route, and the compute route button, so that the web will look more neat.


texts: Loading legend items with hidden label by default in QGIS

texts: Is there a way to set up the loading of the legend items in the layout manager with hidden labels by default? I can't seem to find the settings for that if there are any.


texts: Seeking tutorial to split group of layers in QGIS

texts:
I am trying to split a group of layers, but I can't find anything online about it. I have a layer with a paper clip on it (in the sidebar). I can't figure out how to break it.

Could you link me to a tutorial that would do this?

I don't think I am typing the correct terminology into Google. I am using QGIS


texts: Getting a Python Error with OpenLayers in QGIS 2.18.15

texts:
So I just recently converted to QGIS from pen and paper, and it was all fine until I opened some layers and tried to use OpenLayers for a Google Satellite view. A yellow bar appeared at the top of the screen saying "Python error: An error has occurred while executing Python code: See message log (Python Error) for more details."

Here is the log:

Traceback (most recent call last):
  File "C:/Users/brana/.qgis2/python/plugins\openlayers_plugin\weblayers\weblayer.py", line 100, in addLayer
    self._addLayerCallback(self)
  File "C:/Users/brana/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 182, in addLayer
    self.setMapCrs(coordRefSys)
  File "C:/Users/brana/.qgis2/python/plugins\openlayers_plugin\openlayers_plugin.py", line 231, in setMapCrs
    extMap = coordTrans.transform(extMap, QgsCoordinateTransform.ForwardTransform)
QgsCsException: forward transform of
(110051.440661, 34780.826688)
PROJ.4:  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 +to  +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 +no_defs
Error: latitude or longitude exceeded limits

As you can see, the error code says "latitude or longitude exceeded limits," but even when I zoom in to a reasonable scale (1:2500) I still get the error.

Any fixes?

I am completely new to QGIS and really only know the basic functionality.


texts:
I have a json file downloaded from the ArcGIS website.
You can simply view it in the browser at this address:
https://www.arcgis.com/sharing/content/items/31c5019cf5104960b50f31362d0d8196/data

I would like to load the data into QGIS. It doesn't seem to be compatible...

I used the ogr2ogr command but with no result.


texts: Stop a Cursor loop when only 1 item is selected with Python

texts:
I am trying to create a Cursor that goes through a selection process and stops looping when only one feature is selected. Basically, as the Cursor loops it subsets a selection, and instead of  going through every row of the table I'd like to save time and stop the Cursor when one item is selected. So it would look something like this:

for row in rows:
    ... (selecting process) ...
    if [number of features selected] = 1:
        [stop cursor]
    else:
        rows.updateRow(row)
del row, rows

Any advice? I can't seem to even find any way to stop a cursor other than deleting it but I don't think I could do that type of command from within the cursor itself.


texts: GDAL understanding multiple EPSG codes for a single GeoTIFF

texts:
What does multiple EPSG codes for a single raster image mean ? Is it that the extent spanning across different projection zones.If so how can one generate it using GDAL when an unprojected/GCS image is projected.
C:\Data>gdalinfo MyImage.tif
Driver: GTiff/GeoTIFF
Files: MyImage.tif
Size is 17932, 12384
Coordinate System is:
PROJCRS["Lambert Conformal Conic",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["Lambert Conic Conformal (2SP)",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",42.166666666666,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-89,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",46.6666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",41.3333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (-411299.273570650490001,269487.460171752958558)
Pixel Size = (42.335081965870593,-42.336080785207301)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2020:03:18 07:47:30
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CC (Windows)
  TIFFTAG_XRESOLUTION=300
  TIFFTAG_YRESOLUTION=300
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -411299.274,  269487.460) ( 94d10'42.78"W, 44d28'41.48"N)
Lower Left  ( -411299.274, -254802.564) ( 93d47'42.39"W, 39d45'55.59"N)
Upper Right (  347853.416,  269487.460) ( 84d37' 7.10"W, 44d30'40.70"N)
Lower Right (  347853.416, -254802.564) ( 84d56'35.77"W, 39d47'45.79"N)
Center      (  -31722.929,    7342.448) ( 89d23' 4.29"W, 42d13'55.71"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha


texts: Finding neighbors of latitude longitude pair location

texts:
I have a list of points (more than a few million) defined by latitude and longitude.
I want to find out the neighbors of each point.
For example, I want to know for each point, what the nearest point, or what the nearest 4 points? or how many neighbors are within 10km radius?

I want to do it programmatically, please provide some possible solutions.

So far I have used a for loop to find out the nearest neighbor:

//C++ example
for (i=0;i<npoint;i++)
  {
     lat1=latitude[i];
     lon1=longotide[i];
     for(j=0;j<npoint;j++)
        {
            lat2=latitude[j];
            lon2=longotide[j];
            distance = some_function(lat1, lon1, lat2, lon2);
            vDistance.push_back(distance);
         }
      //now we can find the smallest distance from the container
  }

But this only works for a limited data size.


texts: Apply a polygon filter to CartoJS 4 dataview

texts: Checking this example on CartoJS 4 documentation https://carto.com/developers/carto-js/examples/#example-filter-data-using-drawing-tools, you can apply a circle filter to a source. The problem with this is affects what the user is seeing on the loaded layer. There is a way to apply the same filter but to a dataview? That way the user can draw polygons on top of the layer to get the data without affecting the visible layer. So far the examples only shows how to apply 'bounding box' filters to a dataview.


texts: Creating overlayMaps by querying Carto dataset using geoJson

texts:
I am trying to create an overlayMap by querying a Carto dataset using geoJson.

So far I have:

    var baseMaps = {
        "OSM B&W" : OpenStreetMap_BlackAndWhite,
        "Toner" : Stamen_Toner,
        "Relief Map": Stamen_TopOSMRelief,
        "Satellite Imagery": Esri_WorldImagery
    };

    var geojsonFeature = "http://sspradli.carto.com/api/v2/sql?      format=GeoJSON&q={SELECT * FROM mosquitosamplingsites";

    var bicycleRental = L.geoJson(bicycleRental).bindPopup('Mosquito trap site');

    var overlayMaps={
    "Mosquito Sampling Sites": bicycleRental
    };

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

You can view results here


texts: What does the ArcGIS 3D Analyst parameter 'IsClosed3D' mean?

texts:
I am creating multipatches from existing geometries and finding it hard to work out why IsClosed3D shows some as closed while other apparently closed geometries are shown as not closed.

I have worked out that the ordering of the vertices in different faces has an impact and also that the RingType can affect the result.

If I order all the vertices on the attached geometries to be clockwise, they are both not closed. If I change the order of the vertices one is closed and the other is not closed.


texts: Popup overlays on draw (measure) with OpenLayers 3?

texts:
Using the OpenLayers Measure exemple (http://openlayers.org/en/v3.10.1/examples/measure.html), I'm trying to have only a single measure active at a time.

http://jsfiddle.net/Sinnel/yfdejpre/1/

However, the static measure results appear in a overlay and I have not been able to limit/delete old measures from staying on the map. I have tried different approaches for quite some time, the one giving the "best" results being :

if (map.getOverlays().getArray().length >= 2) {
map.getOverlays().clear()};

However, the result is quite inconsistent...
How should I approach this?


texts: Changing the order of the polygon's nodes using the geometry generator in QGIS

texts: How can I, using the geometry generator, change the order of the polygon's nodes, starting from the most northerly node in QGIS?


texts: Removing pseudo nodes from PostGIS database using ST_ModEdgeHeal()?

texts:
I have used osm2po (5.0.0) and its autogenerated SQL file to load some 25% of Europe.osm.pbf as a network table into my PostGIS database. I have also created an associated nodes table with select pgr_createTopology('network', 0.00001, 'geom_way', 'id');

Since it takes literally days to calculate anything on this dataset with no apparent results, I have simplified it, removing all edges that were of inferior road class (osm2po's attribute clazz > 30), so that my network is largely simplified. This however results in more than 800k nodes (of 1m total in database) that are pseudo nodes (degree = 2), that is, they connect two segments of the same road of the same class and max speed allowed and now are merely unwanted residue from original, more complex topology.

In order to delete those nodes I have been trying to use ST_ModEdgeHeal() as suggested in  How to simplify a routable network? and I intend to use code like the one below.

    select ST_ModEdgeHeal('network', ns.id, nt.id) from pseudonodes p
    left outer join network ns on ns.source=p.id
    left outer join network nt on nt.target=p.id
    where ns.clazz=nt.clazz and ns.kmh = nt.kmh;

It appears, however, that my network is not the network ST_ModEdgeHeal() expects. After I queried SELECT CreateTopology('network', 4326); it created a new blank schema network, whereas I am using schema public where network and nodes tables are stored. I do not know how to combine those two schemas so that ST_ModEdgeHeal() would work. So right now I have an error SQL/MM Spatial exception - non-existent edge 2629192. Understandable, as network schema is filled with empty tables.

Therefore the question is: do I have to do any additional coding, apart from that done automatically with osm2po generated SQL, to set up my network and nodes tables as topology and make ST_ModEdgeHeal() work?

To remove pseudo nodes I have also attempted a solution from this question (Find pseudo nodes in free GIS software), but it did not get any result after 2 days of calculations.


texts: Adding cost from line shapefile using ArcGIS Network Analyst?

texts:
I have a network dataset (a road network), and require a matrix of Cost Distances between many points, which works well.

I have now added a line shapefile (representing rail lines). I wish to use it to add cost when a rail line is crossed. But of course as a line shapefile I am only offered the chance of restriction or scaled cost - but I wish to add a certain value when the rail line is crossed.

I am using Arcmap 10.1 (advanced license).


texts: Openlayers shows and paints polygons in the wrong place on the map

texts:
In the following code I call a geoserver service.

Initialize map

var map = new ol.Map({
        target: 'map',
        layers: [
            new ol.layer.Tile({
                source: new ol.source.OSM()
            })
        ],
        view: new ol.View({
          center: ol.proj.transform([0,0],'EPSG:4326','EPSG:3857'),
          zoom: 2
    })
    });

Then I make a request to get service from geoserver

var sourceVector = new ol.source.Vector({
                // The service returns a GML and also allows return 500 polygons
                format:new ol.format.GML2({}),
                url:'http://localhost:8080/geoserver/siap/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=siap:ppb_curp&maxFeatures=500',
                loader: function (extent, resolution, projection) {
                        //var url = wfsBaseUrl;
                        $.ajax({
                            url: 'http://localhost:8080/geoserver/siap/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=siap:ppb_curp&maxFeatures=500',
                            dataType: "xml",
                            success: function(response) {
                                console.log(sourceVector.getProjection());
                                var gml2format =new ol.format.GML2();
                                sourceVector.addFeatures(gml2format.readFeatures(response,
                                {
                                  featureProjection: 'EPSG:4326',
                                  dataProjection: 'EPSG:3857'
                                }
                                ));
                                console.log(gml2format.readFeatures(response));
                                map.getView().fit(sourceVector.getExtent(), map.getSize());

                                console.log(sourceVector.getExtent());
                            },
                            error: function (jqXHR, textStatus, errorThrown) {
                                alert("The request has failed: " + textStatus);
                            }
                        });
                    }
            });
        //Vector to add styles to the polygon
        var layerVector = new ol.layer.Vector({
                name:'test',
                source: sourceVector,
                style: new ol.style.Style({
                    stroke: new ol.style.Stroke({
                        color: 'rgba(0, 0, 255, 1.0)',
                        width: 2
                    })
                })
            });
        map.addLayer(layerVector);

and all good but the polygons shows them to me by the avid attached image.

The polygons must show them in the area of Mexico, I need some conversion? Any suggestion will help me a lot


texts: Using QGIS C++ API In Custom GIS Application

texts:
I'm trying to develop a custom GIS application using Qt. Which will show a map and have some figure drawing options on the map. While searching on the web, I came across with QGIS. I wonder if I could use QGIS C++ API (QgsMapCanvas etc.) for basic map functionality in my app.
I wrote a simple Qt Application to use QGIS API. It is compiled and linked successfully. But when I tried to run the application it gives: " The program has unexpectedly finished." error. And crashes.
My development environment:

Qt 5.12.0
MSVC 2015
Windows 10 x64
OSGeo4W64 is installed including "qgis-ltr-dev" package  (QGIS version 3.11.0-8) http://download.osgeo.org/osgeo4w/osgeo4w-setup-x86_64.exe (64bit)

My pro file:
#-------------------------------------------------
#
# Project created by QtCreator 2020-06-10T09:28:31
#
#-------------------------------------------------

QT       += core gui xml

greaterThan(QT_MAJOR_VERSION, 4): QT += widgets

TARGET = Hello_QGis_2
TEMPLATE = app

CONFIG += c++11

SOURCES += main.cpp

INCLUDEPATH += C:\OSGeo4W64\apps\qgis-ltr-dev\include

LIBS += -L"C:\OSGeo4W64\apps\qgis-ltr-dev\lib" -lqgis_core -lqgis_gui

DEFINES+=CORE_EXPORT=
DEFINES+=GUI_EXPORT=

# Default rules for deployment.
qnx: target.path = /tmp/$${TARGET}/bin
else: unix:!android: target.path = /opt/$${TARGET}/bin
!isEmpty(target.path): INSTALLS += target

Source Code (main.cpp):
#include <QApplication>

#include <qgsapplication.h>
#include <QtDebug>

int main(int argc, char *argv[])
{
    QgsApplication app(argc, argv, true);
    qDebug() << "*** QGS Application Instance Created *** ";

    qDebug() << "*** setPrefixPath to C:\\OSGeo4W64\\apps\\qgis-ltr-dev *** ";
    QString PREFIX_PATH = "C:\\OSGeo4W64\\apps\\qgis-ltr-dev";

    QgsApplication::setPrefixPath(PREFIX_PATH);

    qDebug() << "*** initQgis ***";
    QgsApplication::initQgis();

    int exitCode= app.exec();

    QgsApplication::exitQgis();

    qDebug() << "*** Finished *** ";

    return exitCode;
}

I also used following bat file to start Qt Creator:
(As mentioned in QT C++ application using QGIS API)
ECHO Setting up QGIS DEV ENV

set PYTHONPATH=

set OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"

@set QMAKESPEC=win64-msvc2015

@set PATH=%OSGEO4W_ROOT%\bin;%OSGEO4W_ROOT%\apps\qgis-ltr-dev\bin;%PATH%

@set INCLUDE=%INCLUDE%;%OSGEO4W_ROOT%\include;%OSGEO4W_ROOT%\apps\qgis-ltr-dev\include

@set LIB=%LIB%;%OSGEO4W_ROOT%\lib;%OSGEO4W_ROOT%\apps\qgis-ltr-dev\lib

path %OSGEO4W_ROOT%\bin;%SYSTEMROOT%\System32;%SYSTEMROOT%;%SYSTEMROOT%\System32\wbem;C:\OSGeo4W64\apps\Qt5\bin;C:\OSGeo4W64\apps\Qt5\Tools\QtCreator\bin;%PATH%

set VS140COMNTOOLS=C:\Program Files (x86)\Microsoft Visual Studio 14.0\Common7\Tools
call "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" x64

start "Qt Creator" /B D:\Development\Qt\Qt5.12.0\Tools\QtCreator\bin\qtcreator.exe %*

As I read from the documentation, the library files for Windows (*.dll) in OsGeo4W64 produced using MSVC 2015. And I tried to use those libraries with the same compiler. So I can not  figure out what am I missing?


texts: Merging several ASCII (.asc) files together and then converting it to a GeoTIFF (.tif)

texts:
In QGIS I'm trying to merge several ASCII (.asc) ESRI files in REGCAN95 (EPSG:4083) projection together to make one big file, and then convert it to GeoTIFF (.tif) in WGS84 (EPSG:4326) projection for use in an orthographic scenery generator (Ortho4XP).
So far I have tried importing the .asc files as raster layers (Layer>Add Layer>Add Raster Layer). I then use (Raster>Miscellaneous>Merge), select all the layers and save it as a .tif. I then use (Raster>Conversion>Translate) with the output projection set to WGS84 (EPSG:4326).
However when I attempt to use the resulting file in the scenery generator, it generates a completely flat landscape as if the elevation data has been lost. I can't figure out what I'm doing wrong.
Does anyone have any suggestions for how to do this in QGIS, or using an alternative method to achieve this?


texts:
I use OSGeo4W on Windows 7 (64 bit) and included libpng16.dll
When I try

C:\OSGeo4W64\bin>gdal_polygonize star.png bla

I get the following error:

ERROR 4: `star.png' does not exist in the file system,
and is not recognised as a supported dataset name.
Unable to open star.png

I copied star.png into C:\OSGeo4W64\bin (this is where gdal_polygonize batch file, python file and compiled python file are located )

What is the mistake?

UPDATE: I found out that I get the same mistake with a TIFF file but it works on JPEGs. Maybe the library/module for PNGs can't be found?


texts: how to pass coordinates to a variable from a geojson source in mapbox gl

texts:
I tried with this code but it wont load the coordinates.the source is fine because it appears in the map when i added it..

map.addSource("villa-route-source",{
type: "geojson",
data: "https://api.myjson.com/bins/h71ul"
})
var route = map.addLayer( {
"id": "villa-route",
"type": "LineString",
"source": "villa-route-source",
"paint": {
       "line-color": "#888",
       "line-width": 8
     }
});
var coordinates = route.features[0].geometry.coordinates;

my purpose is to use the geojson into the turf.along and take its coordinates at specific interval and make a line animation..
I can easily put the geojson code easily and it would work but i want it to be a external source.


texts: ERROR 999999: Error executing function - Failed to execute (CopyRaster)

texts:
ERROR 999999: Error executing function
  Unexpected operation
  No spatial reference exists

This is a fragment from my log report:

  ERROR 999999: Error executing function.
  Database user name and current user schema do not match. [workspace.admin.StreetLightOutages]
  No spatial reference exists.
  DBMS table not found
  Failed to execute (CopyRaster).

I can't import a raster to sde schema using another user (not sde user).
The used user has all privileges (admin privileges) but an error occurred while import  raster dataset to postgresql database using ArcGIS for desktop.
While importing the same raster dataset to the same database using sde user, the operation performed successfully.

Can you help me please whats is wrong in my actions?


texts: Where to find open tabular hyperspectral data for machine learning research?

texts: Where can i find open datasets regarding " Wheat/Maize Leaf Spectra with Associated water content, Potassium and Nitrogen Measurements" for machine learning research i.e. hyperspectral tabular data collected from a field spectroradiometer along with chemical parameters.


texts: Fast data representation on the vector layers

texts:
I have found an interesting article about how HTML5 supports representation on large datasets.
Article can be found here. In the article it is said that HTML5 and feature simplification can increasingly improve visualization speed. Does anybody knows about the source code (or just explain what is needed to be done to reach such results) or suggest how OpenLayers supports feature simplification?


texts: I saved a map, and when I went to open it back up, I was prompted with a dialogue box saying Handle Bad layers.  What does this mean? The map was fine when I saved it.


texts: Consuming ArcGIS Online map service using static IP address for VPN?

texts:
I'm currently using the URL http://services.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer for my iPad map service.

However, we need to connect to VPN to access the internal network at the same time, and the VPN server blocks the public Internet, and only allows me to setup IP Address based public Internet access.

Is there static IP addresses associated to the above URL that I could specify?


texts: Google Maps: obtaining the border of an area

texts:
When using the geocode API it is easy to obtain the (centre of the) coordinates of the address. For example: Rotterdam. If you look up the same location in Google Maps you also get the surface area of Rotterdam.

I would like to do the same using the API, is this possible?

If yes, how?


texts: How to connect SQL SSMS to ArcGIS file geodatabase

texts:
I've found this reference article talking about how to query a file geodatabase using SQL.

http://resources.arcgis.com/en/help/main/10.2/index.html#/Using_SQL_for_reporting_and_analysis_on_file_geodatabases/018s0000000t000000/

My problem is that it doesn't provide any details of how to to connect SSMS (or other query tool) to the file.  Is a driver required? Linked server?

Normally any queries about this result in me directed to ArcSDE. But his article suggests file Geodatabases can be queried by SQL.


texts: definition query and mosaic dataset from gdb ArcMap 10.1 SP1

texts:
I have a Python AddIn Combobox to select a value from the drop-down list. This value changes the definition query of the layers in mxd. This works fine for "normal" feature classes, not for a mosaic dataset. On the Layer Properties dialog box the DefintionQuery tab is empty.

layer_lst = arcpy.mapping.ListLayers(arcpy.mapping.MapDocument("CURRENT"))

for layer in layer_lst:

   if layer.name in ('dted_hillshade','dted_dem6x6') and layer.supports("DEFINITIONQUERY"):

      layer.definitionQuery = "\"NAME\"=\'000E00N.DT2\'"
      arcpy.RefreshActiveView()

Does anyone one know a workaround to set a correct Definition query for mosaic dataset?


texts: How to get the "modified" date for a FeatureClass in a file geodatabase via ArcPy?

texts:
There is a closed question which is similar to or identical to this question: Getting date modified for a feature class within a file geodatabase via arcpy

Looking at a file geodatabase featureclass in ArcCatalog, there is a Modified attribute showing the date that the featureclass was last changed:

This is not the same as the Modified date for its parent file geodatabase (which means that I can't obtain the featureclass 'modified' value via the operating system):

I'm limited to using the older ArcMap/ArcCatalog version of ArcPy, and I can't see where to find the 'modified on' date in the documentation, eg http://desktop.arcgis.com/en/arcmap/latest/analyze/arcpy-functions/featureclass-properties.htm

How can I access the Modified date via ArcPy?


texts: Creating a QGIS Plugin with Plugin Builder - recipe from QGIS Python Programming Cookbook by Joel Lawhead - QGIS 2.18

texts:
I have a trouble with recipe "Creating a QGIS Plugin" from QGIS Python Cookbook by Joel Lawhead. I use QGIS 2.18, Plugin Builder and OSGeo4W Shell.
I can't walk through this step:

  Compile the ui file using the pyuic4 tool:

pyuic4 -o ui_myplugin.py ui_myplugin.ui

OSGeo4W Shell displays an error:

  No such file or directory: "ui_myplugin.ui"

What is the cause?


texts: How to get information about layer-structure on a WMS

texts:
I'm a complete beginner in the world of GIS. And I've been wondering: how do I know where a WMS is, and what layers there are that I can import into my map?

To make it a bit more concrete: I just found the wonderful OpenLayers tutorial here. In this tutorial we import a layer from a WMS by coding:

var vectorLayer = new ol.layer.Tile({
        source: new ol.source.TileWMS({
        preload: Infinity,
        url: 'http://felek.cns.umass.edu:8080/geoserver/wms',
        serverType:'geoserver',
        params:{ 'LAYERS':"Streams:Developed",
            'TILED':true}
    })
});

Now this works fine, and surely is done like that in lots of tutorials.
However, what I couldn't find on any tutorial I've seen, is how do we actually know about the data we are importing? The above example assumes that we know that there is a WMS behind the url http://felek.cns.umass.edu:8080/geoserver/wms, that we are dealing with a GeoServerand not, for instance, with a MapServer, and, most amazingly, that we know that there is a layer called Developed in the Streams workspace on that server.

Surely there is something that I am missing. Where does one get all that information? I couldn't find it on the wesites of WMS-providers like this one.


texts: Can osm way records in SQL spatial be merged to one geographic record if they have same Id?

texts:
I've imported England shp file of roads to mssql using OSM2MSSQL and I've noticed that the same road has couple of records with same Id but different geographic (probably it was segmented).

Is there a way to convert these records to one record with entire geographic data? I want to have unique Id in my streets table instead of having two tables, one with unique Id and one with geographic data.


texts: updating a feature class from an EXCEL spreadsheet

texts:
I created feature classes that only contains the data from an EXCEL file.  I did this to enable the search function to search the data in these tables.  The FCs also have relationships between themselves and other feature classes that have spatial reference.  The database is an ARCSDE SQL Express (Parent) on a server and replicates to three other servers (all Child) using ArcGIS 10.1.

The problem: I receive routine updates to the data in the form of an edited version of the original EXCEL file.  There are over 30,000 records with scores of fields.  The edits can exceed one hundred per day.  I tried running delete rows and upload the new data on each feature class with it's respective spreadsheet, but the delta tables do not pick up those changes so the changes do not replicate to the child.  I have explored Data Management tools for "Compare" and tried running Oject Loader during an edit session - all to no avail.

Any suggestions?


texts:
Oddly, it should be simple, but I can't connect to a local PostGIS database.  My installation is with Mac OS X (10.8), Postgres.app (in doc it's said to have PostGIS included) and Quantum GIS 1.8.

When I want to add a connection to Postgres, I'm specifying name, host (localhost), port (5432), database name, user name and password (service field left empty).

I can connect to a database with Navicat, so database server seems to run correctly.  With the exact information in QGIS used in Navicat, it can't connect... why?

Thank you!


texts: QGIS geometry error:Hole lies outside shell

texts:
I am using QGIS 3.8 and I am using a layer from OSM of Tokyo city.
The layer is consisted from many polygons and I have tried to merge it but I got error.
This is how the original layer was look like:

This is the parts I manage to merge:

When I have checked the geometry validity no error appeard on the problematic areas, beside this one that appeard on area that I manage to merge. I have tried to check but couldn't find any visible error:

I really don't know how to fix it and how to merge the "green layers" with the pink ones. I have seen previous messages in this forum but they were using scripts.

any idea ?


texts: QGIS Expression Dialogue for Reclassify values problem

texts:
I'm trying to reclassify values (simple) for a raster grid within a Model Designer (process modeller). I want to incorporate a 'number input' within the range for reclassifying e.g.:
[0, @MaxDepth, -99999]
Instead of [0, 0.1, -99999]
But this expression does not work and says it is Invalid. It is expecting an array of values with is what I thought [x, y, z] format was.


texts: ArcGIS Polygon to Raster Cell Assignment Manual Threshold

texts: I was wondering if there was any way to apply a user-set area threshold for the cell assignment parameter of the Polygon to Raster tool in ArcGIS. For example, I'd like to only assign a cell if the combined area of the input polygon is greater than 80%. Currently, if you're using "Maximum combined area," the tool just assigns the cell to whichever is the maximum, i.e. greater than 50%, but I was hoping to increase that threshold.


texts: Overlay shapefile data with polygon shapefile

texts:
I would like to add the overlapping data of a worldwide shapefile, to the different European nuts3 regions (polygon shapefile). However, there is something wrong with the script below as it does not work on the shapefiles that I am using.

### download the files available through this link:
https://drive.google.com/open?id=0By9u5m3kxn9yUy1xVDF2NV9TajA

#set your working directory
nuts <- readOGR(".", layer = "NUTS_RG_60M_2010")
grid <- readOGR(".", layer = "a3_mean_annual_runoff_1950_2000")

#projectstring?
#Do I have to do something here to make sure the projections are fine?

#then:
library(raster)
new_spdf <- intersect(grid,nuts) # (this normally works but it does not with these files?

grid_nuts <- gIntersects(new_spdf,nuts,byid = TRUE)
#Take all the values in a particular NUTS polygon, multiply by each the areas of the
#corresponding new grid cells then divide by the total area of the NUTS polygon

#Added the if statement to avoid non-overlapping NUTS polygons

for(i in 1:length(nuts)){
  if(any(grid_nuts[i,])){
    nuts@data$average_spatial_value[i] <- mean(new_spdf@data$value[grid_nuts[i,]]*
                                                 gArea(new_spdf[grid_nuts[i,],])/
                                                 gArea(nuts[i,]))
  } else {
    nuts@data$average_spatial_value[i] <- NA
  }

}

***EDIT : I installed the latest version of R (see sessionInfo() below). The error message I receive is the following:

> new_spdf <- intersect(grid,nuts)
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :
  TopologyException: Input geom 1 is invalid: Self-intersection at or near point 13.3081265 48.608032999999999 at 13.3081265 48.608032999999999
In addition: Warning message:
In intersect(grid, nuts) : non identical CRS

> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgeos_0.3-21 raster_2.5-8 rgdal_1.1-10 sp_1.2-3

loaded via a namespace (and not attached):
[1] tools_3.3.1     Rcpp_0.12.7     grid_3.3.1      lattice_0.20-33


texts: UI CHARTS ( for indices time series in python API of Google earth Engine)

texts: If I use python api for google earth engine and want to make a UI chart for mean value of NDVI for a particular year , how does the UI.CHART function work in python api for google earth engine.


texts: Processing drone Images of body of water?

texts:
I took some pictures of a body of a body of water using a multi-spectral sensor attached to drone. When trying to process the image using Pix4D, more than 50% of pictures were unable to calibrated.

I contacted their customer support to see is there was something that I was doing wrong and their answer was that "images that capture only water, cannot be calibrated (optimize internal and external orientation of the camera). The reason is that it is very hard to find key-points on water surfaces (water is very uniform)".

I am wondering if there is something that I can try in order to stitch all of my images.


texts: Generating SQL expression from multivalue pick list in ArcGIS tool

texts:
I have created a script tool that allows users to select a single value from a pick list, which gets inserted into a query expression.

The pertinent part of my code is pretty simple:

Input = arcpy.GetParameterAsText(0)

arcpy.MakeQueryTable_management(inTable, "TabTemp", "", "", "", "FieldName = " + str(Input))

So if a user selects as input "Item1", then the resulting SQL expression is FieldName = Item1.

I am hoping to make the input a multi-value pick list, so that a user could potentially select one, or multiple values as Input.  For example, if the user selects Item1, Item2 and Item3, then the SQL expression should look like:

FieldName = Item1 OR FieldName = Item2 OR Field Name = Item 3

I'm not quite sure how to get started with the python code to do this.  Can anybody offer suggestions for creating the necessary code to compile the SQL expression from an unknown number of values obtained from a multi-value pick list?


texts: Encode CSV to utf-8 before publishing to Enterprise Portal

texts:
I have a two step process here that's not working. 1) A script successfully writes some requested *.json to a *.csv file. 2) The file is then supposed to be added then published to Enterprise Portal, but is failing at this point with a KeyError. My earlier post shows that once the *.csv is set to utf-8 it publishes without the KeyError. But, I'm lacking a way to do this pythonically.

(run as a stand-alone script)

import json,requests
import os,sys
import csv

path_to_file = os.path.join(r'C:\Users\jpm\Downloads', 'c19_Vaccine_Current.csv')
#idph url goes here
idph_data = 'https://idph.illinois.gov/DPHPublicInformation/api/covidVaccine/getVaccineAdministrationCurrent'

##get json data from idph
response = requests.get(idph_data,verify=True)
#read the json response and keep the VaccineAdministration part
data = response.json()['VaccineAdministration']

#write to file
with open(path_to_file, 'w', newline='', encoding='utf-8') as csvfile:
    f = csv.writer(csvfile)
    #write the headers of the csv file
    f.writerow(['County','AdminCount','AdminCountChange', 'RollAvg', 'AllocDoses', 'FullyVaccinated',                    'FullyVaccinatedChange', 'ReportDate', 'Pop', 'PctVaccinated', 'LHDInventory', 'CommInventory',
                'TotalInventory', 'InventoryDate'])
    for elem in data:
        #get the values for all the keys (i.e. CountyName, AdministeredCount, etc...)
        f.writerow([elem['CountyName'], elem['AdministeredCount'], elem['AdministeredCountChange'],
                    elem['AdministeredCountRollAvg'], elem['AllocatedDoses'],
                    elem['PersonsFullyVaccinated'], elem['PersonsFullyVaccinatedChange'],
                    elem['Report_Date'], elem['Population'], elem['PctVaccinatedPopulation'],
                    elem['LHDReportedInventory'], elem['CommunityReportedInventory'],
                    elem['TotalReportedInventory'], elem['InventoryReportDate']])

(run from a Pro Notebook)

from IPython.display import display
from arcgis.gis import GIS
import os
gis = GIS('Home') #uses current Notebook as workspace

csv_file = r'C:\Users\jpm\Downloads\c19_Vaccine_Current.csv' #path to CSV
csv_item = gis.content.add({}, csv_file) #add CSV to Enterprise Portal
display(csv_item) #display it here

#location parameter needs to be set to "none" when publishing to Enterprise Portal
params={"type":"csv","locationType":"none"}
csv_item.publish(publish_parameters=params) #publish to Enterprise Portal

#results in this error!
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
In  [7]:
Line 3:     csv_item.publish(publish_parameters=params) #publish to Enterprise Portal

File C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\gis\__init__.py, in publish:
Line 10353: elif not buildInitialCache and ret[0]['type'].lower() == 'image service':

KeyError: 'type'
---------------------------------------------------------------------------

If I open the *.csv after the 1st step with this
with open(r'C:\Users\jpm\Downloads\c19_Vaccine_Current.csv') as f:
    print(f)

it prints the encoding as cp1252-- which I've read is because open() uses my machine's environment to choose the appropriate encoding.
I've set the encoding to utf-8 in the original with open() as csvfile. So, how do I bridge the gap? How do I make sure the *.csv is encoded as utf-8 in step 1 before I try publishing it to Enterprise Portal in step 2?


texts:
Wondering why gdal2tiles is throwing errors, I've added the gdaldata and gdalplugins environmental variables and my python is at Python27/ArcGIS10.3/, and am wondering if that is messing with the OSGeo4W python located in the C:\OSGeo4W\bin folder, or do I need to direct my environment variables every time I want to use OSGeo4W vs. ArcGIS python versions?

Or is it working and the error I'm seeing below has to do with something else totally separate?

OSGEO4W home is C:\OSGeo4W\

                   -={ OSGeo4W Shell Commands }=-

       bmp2tiff             listgeo              raw2tiff
       cs2cs                nad2bin              rgb2ycbcr
       curl                 ncdump               rpcgen
       dllupdate            ncgen                rpcinfo
       fax2ps               nearblack            testepsg
       fax2tiff             nircmd               textreplace
       gdaladdo             nircmdc              thumbnail
       gdalbuildvrt         ogdi_import          tiff2bw
       gdaldem              ogdi_info            tiff2pdf
       gdalenhance          ogr2ogr              tiff2ps
       gdalinfo             ogrinfo              tiff2rgba
       gdallocationinfo     ogrlineref           tiffcmp
       gdalmanage           ogrtindex            tiffcp
       gdalserver           opj_compress         tiffcrop
       gdalsrsinfo          opj_decompress       tiffdither
       gdaltindex           opj_dump             tiffdump
       gdaltransform        osgeo4w-setup-work   tiffinfo
       gdalwarp             osgeo4w-setup        tiffmedian
       gdal_contour         pal2rgb              tiffset
       gdal_grid            portmap              tiffsplit
       gdal_rasterize       ppm2tiff             w9xpopen
       gdal_translate       proj                 xmlcatalog
       geod                 psql                 xmllint
       geotifcp             python               xxmklink
       gif2tiff             pythonw
       gltpd                ras2tiff

       epsg_tr           gdalmove          make-bat-for-py
       esri2wkt          gdal_auth         mkgraticule
       gcps2vec          gdal_calc         o-help
       gcps2wld          gdal_edit         o4w_env
       gdal2tiles        gdal_fillnodata   pct2rgb
       gdal2xyz          gdal_merge        rgb2pct
       gdalchksum        gdal_polygonize   setup-test
       gdalcompare       gdal_proximity    setup
       gdalident         gdal_retile
       gdalimport        gdal_sieve

GDAL 1.11.2, released 2015/02/10

C:\Users\maugust\Desktop>gdal2tiles.py C:\Users\maugust\Desktop\example_tifs\Pie
rce_2010_84.tif C:\Users\maugust\Desktop\example_tifs\tilegdal
Traceback (most recent call last):
  File "C:\OSGeo4W\bin\gdal2tiles.py", line 2278, in <module>
    gdal2tiles.process()
  File "C:\OSGeo4W\bin\gdal2tiles.py", line 482, in process
    self.open_input()
  File "C:\OSGeo4W\bin\gdal2tiles.py", line 856, in open_input
    self.out_ds.SetMetadataItem('NODATA_VALUES','%i %i %i' % (self.in_nodata[0],
self.in_nodata[1],self.in_nodata[2]))
IndexError: list index out of range

C:\Users\maugust\Desktop>


texts: arcgis 10.1 account manager, authorizeSoftware script

texts:
I am trying to run "authorizeSoftware" file to complete creating new account manager in Arcgis 10.1, but after running

./authorizeSoftware

I got this message:

  -------------------------------------------------------------------------- Starting the ArcGIS Software Authorization Wizard

  Run this script with --help for additional information.
  -------------------------------------------------------------------------- Unable to access X server (:0.0).  Please check your DISPLAY variable
  or run this script in silent mode.  Run with the --help option for
  more information

Any idea how I can fix this?


texts: ST_Value, multiple rasters - need unique columns for rast values

texts:
I have a grid of cells that I am converting to centroids to 'sample' raster values in the same locations as the centroids. All rasters (4) are in the same table with same resolution and SRID.

I want to return rast values for each surface at each point. The rast values need to be placed in their own unique column based on the value in the 'filename' column and aliased to something more useful than 'rast' - slope, aspect, elev, and tpi. Each centroid will be one record and will have values in each of the aliased rast columns with no empty fields.

This query gets close:

WITH centroid AS (
 SELECT gid, ST_Centroid(grid.geom26912) as centroid_geom
 FROM production.grid
 WHERE field = 'sre'
 AND resol_ft = 60
)

SELECT DISTINCT gid, centroid_geom,
 CASE WHEN filename
  LIKE 'sre10_a%'
  THEN ST_Value(rast, centroid_geom)
 END AS aspect_deg,
 CASE WHEN filename
  LIKE 'sre10_s%'
  THEN ST_Value(rast, centroid_geom)
 END AS slope_deg,
 CASE WHEN filename
  LIKE 'sre10_e%'
  THEN ST_Value(rast, centroid_geom)
 END AS elev_m,
 CASE WHEN filename
  LIKE 'sre10_t%'
  THEN ST_Value(rast, centroid_geom)
 END AS tpi
FROM production.surface10m, centroid
WHERE ST_Intersects(centroid_geom, rast)
 AND filename LIKE 'sre10%'
ORDER BY gid;

However, this results in each gid being repeated 4 times, once for each raster. In each of the four columns there is one (aliased) rast value and 3 NULL values:

My question is how can I avoid the repeating GID and combine the rast values into one line?


texts:
I learned about supervised classification using random forest algorithm.

Details:

Classifier: smileRandomForest with 10 decision trees.
Classes: forest, vegetation non-forest, non-vegetated area, and water.
Class Properties:  name and id.
Sampling: sampleRegions with 500 limit for each class and also, I used rancomCloumn function for each class sample.

Then, I divided my total samples into two classes:

For classification (70% from total samples)
Second, for validation (30% from total samples)

Then, I got my classification result.
After that, I tried to validate it with validation samples (30%)
and I got 100% accuracy.

Is that normal? Or perhaps I made some mistakes in my code.

This my GEE code: https://code.earthengine.google.com/82769aa8927ccc5262eedfa24c657799


texts: How can I create a kml file from a geotiff with alpha layer

texts:
I am trying to create a kml version of a map tile using gdal_translate on a geotiff with alpha layer (link here) as follows:

gdal_translate -of KMLSUPEROVERLAY input.tif output.kml

however, the resulting kml file has a black area where the alpha layer should be - see here for a link to all the input and output files (and subdirectories produced by the operation).

How can I produce a kml file without the black area (also how I could generate only a single standalone .kml file, rather than a .kml linking a bunch of subdirectories as in the zipfile linked)?

gdalinfo on the input geotiff:

Driver: GTiff/GeoTIFF
Files: input.tif
Size is 3822, 3837
Coordinate System is:
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],
        UNIT["degree",0.0174532925199433],
        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"]],
    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"]]
Origin = (-13539823.138956749811769,4690722.589982722885907)
Pixel Size = (0.122995666356360,-0.122995666356360)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DOCUMENTNAME=/mnt/dd/data/1417212380_DDSPIDER/tmp-0-0.tif
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-13539823.139, 4690722.590) (121d37'49.08"W, 38d47' 0.99"N)
Lower Left  (-13539823.139, 4690250.656) (121d37'49.08"W, 38d46'49.09"N)
Upper Right (-13539353.050, 4690722.590) (121d37'33.88"W, 38d47' 0.99"N)
Lower Right (-13539353.050, 4690250.656) (121d37'33.88"W, 38d46'49.09"N)
Center      (-13539588.094, 4690486.623) (121d37'41.48"W, 38d46'55.04"N)
Band 1 Block=3822x1 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=3822x1 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=3822x1 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=3822x1 Type=Byte, ColorInterp=Alpha


texts: Converting from DGN v7 to shapefile by GRASS plugin

texts:
I read two comments in several topics about conversion from DGN files  to Esri shapefile

Both of these comments are interesting for me, but I seek to have them explained in more detail.

I have DGN file v7. I would like to convert it to ESRsri shapefile to work on QGIS with all attributes.

I tried to use the GRASS plugin, but it does not support DGN file to load into this, so how can I do that as someone said in his comment

First comment:

  v7 is the most current .DGN support. Have you tried loading it into QGIS and exporting it from QGIS? You can do this by right clicking on the .DGN file and choosing Save As, but again, only up to .DGN v7

How can I find save as on the file? I did not understand this point.

Second comment:

  I have been wanting to do the same thing, with the same issues. I have successful imported my DGN files into QGIS using the GRASS plugin. I used the v.in.ogr tool in Grass. It allows you to bring DGN files in as lines, points and polys.

Link :Converting from dgn to shp or dwg files using QGIS?

i could not find not v.in.ogr


texts:
I'm correcting a Land Use polygon, interpreting a drone image using QGIS. When I use the cut tool and merge tool to edit, the polygon leave a smalls polygons like a line, or sometimes leave a vertex from other polygon neighbor outside, like a line too. This issue won't let me work fast. Is there a tool or configuration to avoid this?

this is an example of this situation, thats white line is actually a thin Polygon


texts: ArcMap 10.0 Python script tool - Have the script tool dialog box stay open past successful completion of script

texts: Is there a way to have the script dialog box have the radio button that says "Close this dialog when completed successfully" unchecked for any user who runs the script tool?  Is there configuration somewhere to do this?


texts:
I am using update cursor to change the location of a point with FID = 0 every second, but whenever I run the code my ArcGIS freezes and point doesnot move. I am visualizing an output where the point move vertically from (-74.668,-100) to (-74.668,100). Here is my sample code:

import time
import arcpy
#my workspace is set to a folder containing "WorldPoint.shp"
with arcpy.da.UpdateCursor("WorldPoint",("SHAPE@XY",),'"'+'FID"'+' = 0') as cursor:
     for row in cursor:
         for i in range(-100,100,1):
             yvalue = i
             row[0] = (float(-74.668),float(yvalue))
             cursor.updateRow(row)
             arcpy.RefreshActiveView()
             time.sleep(1)

I am new to Python 2.7.


texts:
I'm working on a map with forest trails and I want to get an elevation profile of them all individually.

I have gone through and interpolated the trail layer with a DEM so that I can use the 3D Analyst tool, I select the trail and create a profile graph. Everything looks beautiful, but I was wondering how I would be able to change the axis measurements.

First off I guess how do you tell what the measurements are (i.e. is it in feet, meters, points)?

Second how do I change that (meters to feet)?


texts:
I am new to ArcMap.

I have polygon layer, which represents buildings, and each polygon has number of inhabitants  in Attribute table. Then, i have a layer which shows noise levels in the area, represented also by polygons.
I would like to compute population exposure to those noise levels. What would be the best way to do that? Is it possible? Each building's facade is affected by more noise levels.


texts: Debugging arcpy.SelectLayerByLocation_management() gives ERROR 000368: Invalid input data?

texts:
When i trying to do Select by location i got :
arcgis scripting.
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000368: Invalid input data.
Failed to execute (SelectLayerByLocation)

I am trying to excecute select by location on current selected feature after select by attribute.

env.workspace = "C:\temp"

input = 'C:\\temp\M0TASHT\M0TASHT.gdb\sLine'
input1= 'C:\\temp\M0TASHT\M0TASHT.gdb\sManhole'
arcpy.MakeFeatureLayer_management('sManhole', 'sManhole_lyr')
in_features = 'sManhole_lyr'
where_clause = '"ManholeSubtype" = ('"4515" ')'
arcpy.SelectLayerByAttribute_management(in_features, "NEW_SELECTION", where_clause)
arcpy.SelectLayerByLocation_management("C:\\temp\M0TASHT\M0TASHT.gdb\sLine", "INTERSECT", "sManhole_lyr", "", "NEW_SELECTION")


texts: Remove the Leaflet zoom control on mobile

texts:
Leaflet includes a simple zoom in/out control on the map, unless otherwise specified:

zoomControl: Whether a zoom control is added to the map by default.

This control is unnecessary on a mobile device, where the device's pinch/zoom gestures work better, but the control is included anyway:

Is it possible to detect whether the user is on a mobile device, and remove the control if so (while leaving the control in place for desktop users)?


texts: ArcGIS 10.1 Network Analyst shortest route starting from one specific point

texts: I would like to create the shortest route starting from one specific point in a network out of 12 different points. I know how to get the shortest route by clicking each point after the other to give arcgis the order. but what i wanna do is to just give the starting point and let arcgis calculate the best/shortest route out of all these points.


texts:
5 years ago I exported all our PDF record drawings to TIFF and georeferenced them using Esri Arc Desktop. This added the .tfw file associated with georeferenced images without modifying the original images. This allowed them to be opened as regular images without being rotated, squared off, and filled in with "no data" pixels around the edges. This is important because it allows for better viewing in image software and edits made to the tiff (notes drawings etc..) don't ruin the georeferencing and still allows them to act like a raster data layer in AutoCAD, Esri, QGIS, and other GIS/Spatial Software.
When I try this with PDF's the result is a PDF that is altered in such a way that it is difficult to view as a standard PDF (it is rotated and squared off with a "no data" pixels around the outside).
Is there a way to reference the PDF so that it acts in the same way an image does when associated with a .tfw file? Can it be done using QGIS or other free/open source product?


texts: Aggregating data and displaying counts with GeoPandas


texts: Sampling areas based on water mask using Google Earth Engine

texts:
I am trying to sample a large image and I only want to keep points that fall within a water mask.
I can successfully sample small numbers of points by using a mask, but the algorithm chokes at larger numbers of pixels. This is surprising because I specify the numPixels argument as a large number, but the resultant featurecollection only consists of a small number of points.
Here is my .sample() call with the masking operation:
var water_mask = ee.Image("JRC/GSW1_3/GlobalSurfaceWater")
                    .select(['max_extent'],['is_water'])

var montana = ee.FeatureCollection("TIGER/2018/States")
                   .filter(ee.Filter.eq("NAME","Montana"))

var my_img = img1.addBands(img2).mask(water_mask)

// Call to .sample()
var samples = my_img.sample({'region':montana,
                          'projection':"EPSG:4326",
                          'scale':30,
                          'numPixels':500})
print(samples.size())
// Even though numPixels is specified as 500, the size of 'samples' is only 7

This behavior is synonymous to sampling an unmasked image and filtering the resultant featurecollection by the mask.
For example,
// First add the water_mask bands to the image, so that they can be sampled
var synonymous = my_img.addBands(water_mask)
                 .addBands(ndvi)

                 // Sample masked image - note numPixels
                 .sample({'region':montana,
                          'projection':"EPSG:4326",
                          'scale':30,
                          'numPixels':500})

                 // Filter feature collection by condition
                 .filter(ee.Filter.eq('is_water',1))

print(synonymous.size())
// This results in only 7 points as well

Is there a more efficient way to do this operation?
Here is a link to my code:
https://code.earthengine.google.com/195b11ed5981b695bfe9bbe2da41acac


texts: whats wrong with my exported raster reclassify map in grass?

texts:
I'm trying to export a reclassify map in grass to a Geotiff format to pass it to DinamicaEgo platform.

I have 6 classes in grass but when I export it, the exported map show me 255 classes and they are numbered from 0 to 255. The classes 1 to 6 are my classification but the other not.

It's a problem when I load the raster to DinamicaEGO because it compare two classified images and make a probability matrix. So I only need the 6 classes not 255 classes.

Any one knows how can I make the exported raster only with the 6 classes?

I used the command r.reclass for reclassify the image and r.out.gdal for export it.


texts:
I have several tables in SQL server. I am trying to visualize my data in ArcGIS. However, I have not been able to do so.

I started off my making a database connection in Arc Catalog to my SQL. I did this correctly and I am able to view my data.

Then; I went to ArcGIS 10.3 -> Make Query layer and selected my SQL connection. I then double clicked on the table I wanted. Checked advanced options and validated and clicked next. I then selected the primary key and tried to select the CRS. I could not because it was greyed out. So I decided to continue and finalized. In the end I am only able to see my connection in ArcGIS but no mapped data.

How can I map my data?


texts: Auto-populate "from" and "to" attribute fields of line layer with Unique IDs of points from point layer?

texts:
I have a line layer where each line segment is snapped to 2 of the points from our point layer in QGIS.  Each point has a Unique ID which I need to carry over automatically to a "from_point" field and a "to_point" field in the line layer attribute table.

Would there be a way to populate the fields after digitizing all the lines and points?

I figure we will need to run some sort of spatial query to relate the 2 points to each line segment and then run some sort of expression to populate the fields.


texts:
I am trying to use this lightning data set. https://www.ncei.noaa.gov/thredds/dodsC/lightning/countPerYear_Aggregation.ncml.html
i do not understand the coordiante system which was used in generating the dataset. I want to understand how to convert a simple lat, long in WGS84 into the coordinate system used to generate the dataset.

PROJCS["NAD83 / Texas Centric Albers Equal Area",

GEOGCS
["NAD83",DATUM
  ["North_American_Datum_1983",SPHEROID
  ["GRS1980",6378137,298.257222101],
  TOWGS84[0,0,0]
],

PRIMEM
["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]],

UNIT["metre",1,AUTHORITY["EPSG","9001"]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["standard_parallel_1",27.5],
PARAMETER["standard_parallel_2",35],
PARAMETER["latitude_of_center",18],
PARAMETER["longitude_of_center",-100],
PARAMETER["false_easting",1500000],
PARAMETER["false_northing",6000000],
AUTHORITY["EPSG","3083"],
AXIS["X",EAST],AXIS["Y",NORTH]]

institution: National Climatic Data Center (http://www.ncdc.noaa.gov)
source: generated using SwdiGridder (version trunk) by Mark Phillips  from input file: ../csv/nldn-2011.csv
history: generated on Fri Dec 28 21:00:49 EST 2012

How can I understand this coordinate system of this dataset?