texts: How to calculate the difference of rasters within the limits where the pixel coordinates are the same?

texts:
I'm writing my raster calculator like Qgis.I'm stuck in one place.
I feed two images to the input, after which I find the boundary of the intersection of the coordinates of the pixels, since my two images differ in coordinates

F=[]
F.append(r'C:\Users\user2\Desktop\1\rastrcalc\pyopencl\2012_2.tif')
F.append(r'C:\Users\user2\Desktop\1\rastrcalc\pyopencl\2011_cut_2.tif')
# For Image Processing

from math import *
import numpy
from osgeo import gdalnumeric
from osgeo import gdal
from osgeo.gdal_array import *
from osgeo.gdalconst import *
import pyopencl as cl
# Set our output file format driver

driver = gdal.GetDriverByName( 'GTiff')
# Set our output files Projection parameters
# from input file number 1

tmp = gdal.Open(F[0])
geoT = tmp.GetGeoTransform()
r1 = [geoT[0], geoT[3], geoT[0] + (geoT[1] * tmp.RasterXSize), geoT[3] + (geoT[5] * tmp.RasterYSize)]
proJ = tmp.GetProjection()

tmp2 = gdal.Open(F[1])
geoT2 = tmp2.GetGeoTransform()
r2 = [geoT2[0], geoT2[3], geoT2[0] + (geoT2[1] * tmp2.RasterXSize), geoT2[3] + (geoT2[5] * tmp2.RasterYSize)]
proJ2 = tmp2.GetProjection()
intersection = [max(r1[0], r2[0]), min(r1[1], r2[1]), min(r1[2], r2[2]), max(r1[3], r2[3])]
>>> print intersection
[286209.3065796, 5180967.071686, 286223.8630979998, 5180947.93927]
b3 = LoadFile(F[0]).astype(numpy.float32)
b4 = LoadFile(F[1]).astype(numpy.float32)

#Reshape the data from 2-D to 1-D

b3r = b3.reshape(b3.shape[0]*b3.shape[1])
b4r = b4.reshape(b4.shape[0]*b4.shape[1])

With the help of variable "intersection" created a frame

Now the question is how to calculate the difference of rasters exactly within the "intersection"
At the location where the pixel coordinates are the same

The principle of the raster calculator Qgis

It is worth noting that finding the frame is taken from here


texts: Changing representation overrides to attribute driven symbology for migration to ArcGIS Pro

texts:
We currently use a lot of symbology representations in our map production workflows in ArcMap and are considering moving to ArcGIS Pro, which does not support representations. Many of them will be ok to make the change using attribute-driven symbology, but I am wondering if anyone else has gone through this and has any suggestions about how to deal with representation overrides.

The override is stored as a blob, but we noticed that we can still query for null in the blob field, so I assume I can use this to pull out every feature that does in fact have an override. But then I don't have much information, if any, about what exactly was changed on each representation symbol. Is there a way to get this information out of the blob field, or some other way to extract this information from the database?


texts:
I have a .bmp image with a map. For creating image i used a Ingeo 4.4 and function MakeImage() What i know:

Height an Width of bmp image
dpi
Map Scale
Image Center's coordinates in meters.

What i want:

How can i calculate some points of image (for example corners) in meters.
Or how can i change a pixel distanse to meters?

What i do before:

For sure i know image center coordinates in pixels:

CenterXpix = Widht/2;
CenterYpix = Height/2;

But what i gonna do to find another corners coordinates. Don't think that:

metersDistance = pixelDistance*Scale;

is a correct equation.

Any advises?


texts: Points to Path using order from array of points


texts: How to display ESRI shapefile in Android

texts:
It's the first time I work with Android and ESRI shapefiles. I tried to find but could not find any tutorials to display ESRI shapefile in Android.

Somebody can help me about this?


texts: Unlist into a single Object - Google Earth Engine

texts:
I'm new to GEE and Java scripts. I'm trying to do something that seems simple, but cannot find any way to do it.

I’ve an image (which is not that big, say 100 km * 100 km) with ~ 12 bands (and a 30 m spatial resolution in this case). My goal is to obtain the mean of each band (in order to center the bands, before a PCA). This is much harder than anticipated, because of user memory limit.

I tried several strategies, from the most naïve one (considering all the pixels and a ‘reduceRegion’ function) to a less computer-intensive one (considering a random subsample of pixels), and now trying to loop it band by band… in order to avoid reaching the user memory limit… while trying to obtain the output format of the initial “reduceRegion” function strategy.

Here is an example : [link] (https://code.earthengine.google.com/11e0e4b9f7ca71c42774f390ddecf0af)

In my last attempt, I've defined a simple function that returns the mean of a given band in a satellite image, and applied it with map() on the multi-band image:

var listofbands = ee.List(['B1', 'B2']);  // say i focus on 2 bands, named B1 and B2

var map_m = function(i) {

  i = ee.String(i)

  var Band_i = image.select(i)  // select the band i in the image

  var Mean_i = Band_i.reduceRegion({ // compute the mean over "sample", with sample e.g. 10 sample points

  reducer: ee.Reducer.mean(),

  geometry: sample.geometry(),

});

return Mean_i} // return this mean

var listOfResults = listofbands.map(map_m).flatten()  // apply the function to listofbands

print(listOfResults)

It returns a "List", which contains one object per band in the image (and each object as the band name and the mean value). So it's not too bad but I do not want a "List", i want a single "Object" with all the band names and mean values in it:

I'm therefore trying to "unlist" this "List" (so that i obtain a single "object"), or convert the list to a single object (Dictionary)... without luck so far. This should be quite straightfoward... any idea?

Best


texts: How to display information when a certain dataset is loaded into QGIS?

texts: I need to show the quality information of a dataset that was manually created by me, that is saved as a text file. The text file will "pop up" when the dataset is added into QGIS. I have no experience in creating plugins in QGIS. Can you help me with the approach to create a plugin for the above mentioned purpose.


texts: How to specify geoid grid for `osr.SpatialReference` instance

texts:
How can I specify e.g. Geoid 12A or Geoid 18 when (for example) converting NAVD88 orthometric heights to ellipsoid heights using the OSR Python bindings?
Of course I would like to be able to specify any .gtx file that is available to OSR.
A working example would be very much appreciated as I'm still a little unsure of how to make this work. I guess it looks something like this:
from osgeo import osr

# Take NAD83/UTM 15N as an example
src_horz = osr.SpatialReference()
src_horz.ImportFromEPSG(4269)
src_horz.SetUTM(15, True)
# Use NAVD88 vertical datum
src_vert = osr.SpatialReference()
src_horz.ImportFromEPSG(5103)
# Define compound SRS
src_srs = osr.SpatialReference()
src_srs.SetCompoundCS('', src_horz, src_vert)

# Define output SRS
dst_horz = osr.SpatialReference()
dst_horz.ImportFromEPSG(4326)
# Use NAVD88 vertical datum
dst_vert = osr.SpatialReference()
dst_vert.ImportFromEPSG(4979)
# dst_vert.SetVerticalShiftGrid('Geoid 12A') ### How do I do this?
# Define compound SRS
dst_srs = osr.SpatialReference()
dst_srs.SetCompoundCS('', dst_horz, dst_vert)

transform = osr.CoordinateTransformation(src_srs, dst_srs)


texts: Subset the shapefile from Sentinel 2 data in Google Earth Engine

texts: I am new to Google Earth Engine and I am working with Sentinel-2 data for salinity mapping. I have the shapefile of my Area of interest, I want that to be subsetted from the data. What commands should I need to get this?


texts: Determining if two geometric figures are approximately equal using PostGIS?

texts:
I would like to know how to determine if two geometric figures (the_geom) are approximately equals.

By the way, I'm working with Python, PostgreSQL and Postgis

I have the following geometries:

0101000020110F0000FE546AD2413E5EC16D56617C93C74FC1

0101000020110F00009D7E6384413E5EC10185DD7F93C74FC1

When I see both figures on the map, the variation is imperceptible.

Any idea?


texts:
How do I get a list of layers by clicking on overlapping WMS layers in the Leaflet map?  The best solution is to have the list in a single popup on the map.

<!DOCTYPE html>
<html>
<head>
    <meta charset="UTF-8">
    <title>WMS Leaflet</title>
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.2.0/dist/leaflet.css" />
    <style>

     body {
            padding: 0;
            margin: 0;
        }
        html, body, #map {
            height: 100%;
            width: 100%;
        }

    .leaflet-popup-content-wrapper {
    width: 420px;
    height: 400px;
    }

    iframe {
    width: 400px;
    height: 350px;
    }
    </style>
</head>
    <body>
    <script src="https://unpkg.com/leaflet@1.2.0/dist/leaflet.js"></script>
    <script src="https://cdn.jsdelivr.net/npm/leaflet.wms@0.2.0/dist/leaflet.wms.min.js"></script>

<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4js/2.5.0/proj4.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/proj4leaflet/1.0.2/proj4leaflet.js"></script>

    <div id ="map"> </div>
    <script>

        var map =
  L.map('map', {
    center: [43.731739, 10.401401],
    zoom: 2,
    // minZoom:12,
  });

var options = {'transparent': true};
var source = L.WMS.source("http://geoportale.regione.lazio.it/geoserver_nsit/ptpr/wms", options);
var layer1 = source.getLayer('zone_umide');
var layer2 = source.getLayer('aree_protette');
var control = L.control.layers({}, {
    'Layer 1': layer1,
    'Layer 2': layer2
})
control.addTo(map);

    </script>
    </body>
</html>


texts: Is it possible to create a standalone model in ArcMap Modelbuilder?

texts: Is it possible to create an algorithm in ArcMap Modelbuilder that could then become a fully standalone applet? I would like to be able to provide this processing algorithm even to users who do not own an ArcMap license.


texts: Reading NV chart data (.bsb and .eap) with QGIS?

texts:
I have some data from NV charts. As it looks like the data is contained in .bsb and .eap files. Is there a possibility to implement these data in QGIS as raster as a background? Loading the .bsb data directly seems not possible. Loading the .eap data works but there is only black and white lines and nothings looks like a map.

This nautical chart data seems to be different from NOAA data as no .kap files exist.

Data can be read using the following software:
nv charts App, Chart Navigator professional, or Chart Navigator Standard.


texts: Data Reviewer Regular Expression for HHMM time format

texts:
I am using Data Reviewer to check to see if a string field is in the time format HHMM using the regular expression check. I think I have the correct logic with this:

^([0[0-9]|1[0-9]|2[0-4])[0-5][0-9]

but data reviewer doesn't pass any valid times with this expression (1234, 1259, 1340, 0013, 0240, 2349 all show up as errors, for example).

Is the Data Reviewer regular expression in SQL format, or is there something I am missing when writing regular expressions for data reviewer? I got the date check in YYYYMMDD format working no problem with:

((17)|(18)|(19)|(20))\d\d((0[1-9])|(1[012]))((0[1-9])|([12][0-9])|(3[01]))


texts: How do I turn a fog belt zone raster map into a geographic map of expected fog?

texts:
I'd like to make a map that shows the amount of fog by San Francisco neighborhood.

I'm using QGIS. I have some of the pieces:

a geoJSON file for each of the neighborhoods.
a raster map with a shapefile for the fog belt zones
a vector file for the contours generated by the fog belts
a map that shows the average number of hours for the Bay Area

My goal would be something like this, where each neighborhood has one value for the average summertime hours of fog & low clouds per day in that neighborhood -

https://i.stack.imgur.com/tKL2x.jpg

I'm not sure how to read the raster map of fog belt zones. My understanding is that the area between each of the contours indicates how much fog is expected in that geographic space.

It's unclear to me how to turn that into something more readily understandable, like the map I linked above.

How do I translate the fog belt zones file into expected fog in geographic space?


texts: how to load QGIS workspace in Qfield. There is no /android folder

texts:
I am trying to load a Workspace into Qfield on my S8 android phone. I have followed the steps on any website i have come across but I am unable to find or access the folder in which my layers and workspace have been placed. When trying to save the folder on my phones internal memory I have put in in the /android/data/ch.opengis.qfield/files/newworkspace but when i try to access this through qfield I cannot seem to get into the /android folder as it is not there on my device. I have pressed the upwards arrow as many times as possible and the first 4 options in my list are as follows

acct

cache

config

cpefs

I have also tried saving the workspace and layers on my sdcard as I can access this but when i get as far as the sdcard/ch.opengis.qfield/files I cannot see the folder in there in which I have created that is holding my workspace information.


texts: Aggregating raster from Combine Tool in Spatial Analyst--->Local?

texts: I have a raster file created through the Combine tool (Spatial Analyst--->Local). The attribute table is showed below. Its resolution is roughly 10m. I would like to aggregate it to 50m with each aggregated cell (including 25 small cells) having the sum of those attributes (SFDVal, MFDVal, MHVal and etc.) from the 25 small cells. I tried several ways (aggregation, zonal etc) but it seems like all of them are working on the actual cell value of the combined raster, instead of other attributes. But I only want to get a sum aggregation on those attributes, not the cell value.


texts: Join spatial point data to polygons in R

texts:
I am trying to perform a spatial join between point data and polygon data.

I have data that indicate the spatial coordinates of an event in my csv file A and have another file, shapefile B, that contains the boundaries of an area as polygons.

head(A)
  month   longitude latitude lsoa_code                   crime_type
1 2014-09 -1.550626 53.59740 E01007359        Anti-social behaviour
2 2014-09 -1.550626 53.59740 E01007359                 Public order
3 2014-09 -1.865236 53.93678 E01010646        Anti-social behaviour

head(B@data)
  code      name                                  altname
0 E05004934 Longfield, New Barn and Southfleet    <NA>
1 E05000448                   Lewisham Central    <NA>
2 E05003149                            Hawcoat    <NA>

I want to join the crime data A to my shapefile B to map the crime events that happen in my area A. Unfortunately I cannot perform an attribute join based code as the code in A refers to different units than the code in B.

I've read a number of tutorials and posts but could not find an answer. I tried:

joined = over(A, B)

and overlay, but did not accomplish what I wanted.

Is there a way to do this join directly or would an intermediate transformation from A to another format be needed?

Conceptually I want to select those points of A that fall into the code areas of B (similar to "join based on spatial location in ArcGIS").

Did someone have this issue and solved it?


texts:
When I try to Import Shapefile into PostGIS using DB Manger Plugin in QGIS the error message show

"Error 3
Creation of data source "TEST_EST_MILION" failed:
ERROR:  function addgeometrycolumn(unknown, unknown, unknown, integer, unknown, integer) does not exist
LINE 1: SELECT AddGeometryColumn('','TEST_EST_MILION',NULL,32640,'MU...
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
"MU...
           ^


texts: OSM labels look as red point in geoserver preview layer. Why?

texts: I've imported shp directory files of OSM to geoserver and the city label layer look like red dots instead of showing labels of street names. Is there a fix for this ? (should I declare this layer kml instead of openlayers ?)


texts: How can I get points from polygon considering a buffer zone where within that buffer only one point will be generated. I have a few polygons for suitable locations. Now I want to make points from them, but I also need to consider at least 500m distance from one point to another. I have tried to densify and vertices to points, generate points along line but none of these are giving expected outcome.


texts:
I want to add HeatMap customized data based in a vector (longitude, latitude, temperature). Currently I have this code:

 data = new ol.source.Vector();
 map = new ol.Map({
     target: 'map_canvas',
     layers: [
       new ol.layer.Tile({
           source: new ol.source.MapQuest({ layer: 'sat' })
       })
     ],
     view: new ol.View({
         center: ol.proj.transform([37.41, 8.82], 'EPSG:4326', 'EPSG:3857'),
         zoom: 4
     })
 });

 // created for owl range of data
 var coord = [message.latitude, message.longitude];
 var lonLat = new ol.geom.Point(coord);
 var pointFeature = new ol.Feature({
     labelPoint: lonLat,
     weight: 20 // e.g. temperature
 });
 data.addFeature(pointFeature);
// create the layer
 heatMapLayer = new ol.layer.Heatmap({
     source: data,
     radius: 50
 });
 // add to the map
 map.addLayer(heatMapLayer);

After all these steps I cannot see any changes in the map.

Can someone help me this issue?


texts: Specify position of curved label in QGIS?

texts:
From this answer in 2012 underdark stated:

  BUT you won't be able to have curved labels anymore if you move/rotate the labels manually, as far as I know.

It's now 2016, and I noticed that I can enable the Curved radio button and specify coordinates for my label (see below). However when I specify these coordinates, my label is no longer curved.

I want to specify a particular orientation for curved label around a circular buffer (which I've converted to a polyline), is there a way to do this?


texts:
I am attempting to write many ASPRS LAS point files into an ESRI File Geodatabase using my compilation of GDAL 1.9.2. The FileGDB driver for GDAL/OGR appears to be incredibly slow when writing large files, taking as much as 45 minutes to write just 8 million point records. Write speeds of FileGDB using GDAL on a SATA3 drive are somewhere on the order of 200 kilobytes per second, which is unacceptably slow when I'm trying to convert terrabytes of data.

I noticed in the FileGDB documentation that defining the FGDB_BULK_LOAD macro should improve performance for large datasets, but I noticed no change in performance when I wrote a line into the "nmake.opt" file with the text "FGDB_BULK_LOAD = YES" immediately after the FGDB_LIB line.

Admittedly, a FileGDB is not the ideal way to store billions of point data records, but that's a gripe for another time. Have I correctly used the FGDB_BULK_LOAD feature? Is that supposed to be in my source code, not the GDAL build?

Thanks.

UPDATE: Proper Usage: (Answered in chat)

The FGDB_BULK_LOAD setting is properly stored as an environment variable for the GDAL/OGR process. This is set on the command line during the ogr exe call as shown by Ragi. Using the GDAL functions, it can be set programmatically for the entire process with

 CPLSetConfigOption("FGDB_BULK_LOAD", "YES");

or just for the current thread using

 CPLSetThreadLocalConfigOption("FGDB_BULK_LOAD", "YES");

FGDB_BULK_LOAD must be set prior to calling FGdbDataSource::CreateLayer(). It was unclear whether OGRCleanupAll() unset this variable, but it is safe to call multiple times to be sure.

Using that option boosted performance to be around 5.5x faster for writing millions to tens of millions of points.


texts: QGIS 2.12.3-1: startup aborts after splash screen "Setting up the GUI"

texts:
Fresh install of QGIS 2.12.3-1 through OsGeo4W64 on a Win7-x64 machine:

When I try to start QGIS through the start menu, through OsGeo4w/bin/qgis.bat or the OsGeo4W Shell with the command qgis I get the same behaviour:

The splash screen appears and just after "Setting up GUI" it closes. No error message appears. Besides QGIS (and its automatic dependencies) no other package was installed through the OsGeo4W setup tool.

Any ideas if it is missing a dependency or how I could troubleshoot this issue with no present error message?


texts: How can I get clean road attributes from OSM data?

texts:
I would like to use osm road data for a project where I need to broadly distinguish between primary, secondary, and tertiary roads.

OSM data seems just fine to work with Arcgis' Network Analyst.

When I download data for WeoGeoMarket, selecting only highways, I face a lot of different classes, sometimes mis-coded (such as 'teriary' instead of 'tertiary').
However, if I visually check on OSM interactive map, I can see that the legend is not messy at all.

My question is therefore the following: is it possible to get OSM road data with clean feature attributes?

Thanks for your help!


texts:
I try to use this code to fill attribute of map with my form but it dose not work?

Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click
        Dim mxdoc As IMxDocument = My.ArcMap.Document
        Dim pmap As IMap
        pmap = mxdoc.FocusMap
        Dim penumlayer As IEnumLayer = pmap.Layers
        If pmap.LayerCount > 0 Then
          Dim player As ILayer = penumlayer.Next
          Do Until player Is Nothing
                Dim pflayer As IFeatureLayer
                If TypeOf player Is IFeatureLayer Then
                    pflayer = player
                    Dim pfclass As IFeatureClass = pflayer.FeatureClass
                    pfclass.Value(pfclass.Fields.FindField("TDate")) = TextBox1.Text
                    pfeature.Value(pfeature.Fields.FindField("Proj_Type")) = ComboBox1.Text
                    pfeature.Value(pfeature.Fields.FindField("Utility_Ty")) = ComboBox2.Text
                    pfeature.Value(pfeature.Fields.FindField("Destrict")) = ComboBox3.Text
                    pfeature.Value(pfeature.Fields.FindField("Team_No")) = ComboBox4.Text
                End If
                player = penumlayer.Next
            Loop
        Else
            MsgBox("NO LAYERS")
        End If
 End Sub


texts: Export ArcMap Symbology to table, text or other format?

texts: In ArcMap how can I export a symbology style to table or text format such as .csv?


texts: Can I adjust the position of individual qgis labels?

texts:
I have been playing with the basic label placement settings, but I'm wondering if there's a way for me to adjust the placement of individual labels?

There are a few that overlap and I'd like to be able to manually adjust them.

Or ... do I have to export an SVG and do that last step there?


texts: Setting Definition Query on ArcPy layer from shapefile for ArcMap

texts:
I am trying to determine options to set a layer definition query using ArcPy.

I know that it is possible when using arcpy.mapping.ListLayers().

However, in this scenario I am loading in shapefiles, and using arcpy.MakeFeatureLayer_management() in ArcPy.

How do I set a definition query on the layer that I create from using this geoprocessor?

I am using ArcMap 10.2.


texts: Case Insensitive WFS Attribute Search with ArcGIS Server?

texts:
I'm using ArcGIS 9.3.1 to host WMS and WFS services.

I have a map (.mxd file) with a layer whose attribute table includes fields with mixed-case string data. When I open this map in ArcMap and perform a Find (Edit -> Find), selecting the appropriate layer and field, I get all matching results back, no matter what case I use in entering the search criteria. Searching for 'fireside', 'Fireside', and 'FIRESIDE' in STR_NAME all return features whose STR_NAME is 'Fireside Ln'.

However, if I take the same map and publish it to WFS in ArcCatalog, when I perform a WFS request against it, I only get matching results when searching for 'Fireside' (not 'fireside' or 'FIRESIDE').

I am performing a WFS request in the form of an XML POST. This is an example request that returns no results:

<wfs:GetFeature service="WFS" version="1.1.0" maxFeatures="1000" resultType='results' xmlns:wfs="http://www.opengis.net/wfs" xmlns:ogc="http://www.opengis.net/ogc" xmlns:gml="http://www.opengis.net/gml">
 <wfs:Query srsName="EPSG:4326" typeName="centerlines:street_centerline">
  <ogc:PropertyName>STR_NAME</ogc:PropertyName>
  <ogc:PropertyName>TYPE</ogc:PropertyName>
  <ogc:PropertyName>Shape</ogc:PropertyName>
  <ogc:Filter>
   <ogc:PropertyIsLike wildCard="*" singleChar="." escape="!">
   <ogc:PropertyName>STR_NAME</ogc:PropertyName>
   <ogc:Literal>*fireside*</ogc:Literal>
  </ogc:PropertyIsLike>
  </ogc:Filter>
 </wfs:Query>
</wfs:GetFeature>

I have also tried adding a matchCase="false" attribute to the ogc:PropertyIsLike tag (based on a discussion I found on osgeo.org), but this also returned no results.

Is there any way to perform a case-insensitive attribute search against a WFS service hosted on ArcGIS Server, either through modifying the request, or by performing some kind of configuration of the layer, map, or service on the server end?


texts: Choosing projection for points if they extend outside of bounds

texts:
I have a project where there are GPS points that are centered around the Chesapeake bay, with some points extending between North Carolina and southern Canada

I have been using the NAD83 Virginia South state plane because I was told that a projected coordinate system is more accurate for distance calculations, but I am unsure if this projection would be reliable with the spread of my data points.

I read that outside of a specific state zone errors increase quickly, but I couldn’t find any information on how error accumulates if points extend outside of the bounds of the state zone (for example, for every 5 km outside of the state zone you’d see error increase by 1 m).

Do you know if this information is available anywhere, or am I misunderstanding how the error works? Would a more general projection for all of North America such as the North America Equidistant Conic be a better choice?


texts: Using arcpy.makequerytable_management with where clause?

texts:
I have a script that works great with selecting dates greater than "my From date" but I can't get the code to work when I want to add a less than to the sql where clause (a to date).

I had add the {0}'".format to get the original where clause to work.  But can't figure out how to add the second date.  I have tried AND and BETWEEN with no luck.

Any idea how to do this?

I'm stuck!

This is the code that works perfectly for the single from date:

where_time_period = "dtrepor > '{0}'".format(object_FromDate)

table_out =  os.path.join(output_gdb, "MyRecords")
if arcpy.Exists(table_out):
    arcpy.Delete_management(table_out)
arcpy.MakeQueryTable_management(incident_table, table_out, where_clause =
where_time_period)

update:
If I run the code with just the where_time begin it works perfectly. If I try to add where_time_end it gives me all the records. I just tried to change the where_clause to this

where_time_begin = "dtrepor > '{0}'".format(object_FromDate) where_time_end =
   "dtrepor < '{0}'".format(object_ToDate)
   arcpy.MakeQueryTable_management(incident_table, table_out, where_clause =
   where_time_begin and where_time_end)

instead of giving me less record it

It returns all the records –

If I try where_clause with "AND" or Between, my code won't run?
   where_time_begin + "AND" +  where_time_end

Here's a bit more of the code

#define workspace#
arcpy.env.workspace = r"C:\Users\mandrews\Documents\ArcGIS\Default.gdb"
dbo_path=r"C:\Users\mandrews\AppData\Roaming\ESRI\Desktop10.5\ArcCatalog\Spillman_rep_OLE_D.odc"
sde_path = "C:\Users\mandrews\AppData\Roaming\ESRI\Desktop10.5\ArcCatalog\RepSpillman.sde"
#mxd_all_incidents = arcpy.mapping.MapDocument(r"C:\AlexaWetmore\monthlyAll.mxd")
mxd_all_incidents = arcpy.mapping.MapDocument(r"C:\Users\mandrews\Documents\monthlyIncidentMaps\monthlyAll.mxd")
output_dir = (r"C:\Users\mandrews\Documents\monthlyIncidentMaps")
output_gdb = (r"C:\Users\mandrews\Documents\monthlyIncidentMaps\monthlyIncidents.gdb")
#call addr temp
addr_temp = os.path.join(r"C:\Users\mandrews\Documents\ArcGIS\Default.gdb", "AdderTemp")
#copy lwmain table with this sql query for from date
incident_table = os.path.join(dbo_path,"dbo.lwmain")
offense_table = os.path.join(dbo_path,"dbo.lwoffs")
where_time_begin = "dtrepor > '{0}'".format(object_FromDate)
where_time_end = "dtrepor < '{0}'".format(object_ToDate)
table_out =  os.path.join(output_gdb, "MyRecords")
if arcpy.Exists(table_out):
    arcpy.Delete_management(table_out)
arcpy.MakeQueryTable_management(incident_table, table_out, where_clause = where_time_begin + "AND" + where_time_end)
#show count of records selected
results = arcpy.GetCount_management(table_out)
count = int(results.getOutput(0))
print (count)


texts:
I'm downloading a point-dataset as an x,y,depth .csv-file. I'd like to interpolate these and use them as an nautical navigational chart.
The problem arises on interpolating the points into a raster (image below). The dataset also contains islands and parts ashore. I'd like to have a method which sort of clips the raster if no points are in the neighborhood. Even better would be a method which doesn't interpolate at these areas.

In the figure below, I indicated the part in green which yields the perfect result: no interpolation and thus no_data. The parts in red indicate the areas which I want to get rid of. Is there a good method or workflow to yield this desird result?

Edit:
Here is a small sample dataset without islands. I'm still preparing a set with islands in it. Samples


texts: How can I plot scatter plots over a shapefile using Python?

texts:
This is my code:

import matplotlib.pyplot as plt
import json
import time
import geopandas as gpd
from descartes import PolygonPatch
from shapely.geometry import LineString

borough = gpd.read_file('London_Borough_Excluding_MHW.shp')
borough = borough.to_crs({'init': 'epsg:4326'})

fig = plt.figure(figsize=(20, 8))
ax = borough.plot()

x, y = ax(df2['lon'].values, df2['lat'].values)
ax.scatter(x,y, marker="*", color='r', alpha=0.7, zorder=5, s=2)
plt.show

This is the error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-40-5f23b6559223> in <module>
     12 ax = borough.plot()
     13
---> 14 x, y = ax(df2['lon'].values, df2['lat'].values)
     15 ax.scatter(x,y, marker="*", color='r', alpha=0.7, zorder=5, s=2)
     16 plt.show

TypeError: 'AxesSubplot' object is not callable

and this is the result I got:


texts: Importing GEE feature collection output as sf object in R

texts:
When you export a feature collection (as a CSV) from Google Earth Engine, the geometry information is stored in WKT format under the variable .geo. The general format for each cell in the vector .geo looks something like this:
{""type"":""Point"",""coordinates"":[-16.129853523180806,12.1466390262208]}
If I read in the CSV of the feature collection:
centroidData1 = fread("Tiles/Demo/count_1_plus_rasterdta.csv")
theoretically, we should be able to use .geo to transform centroidData1 into an sf object using the st_as_sf function.
myCRS= "+proj=longlat + datum=WGS84"

centroidPoints1= st_as_sf(centroidData1,
                          wkt=".geo",
                          na.fail=T,
                          crs=coordRS)

But when I execute this code, R returns the following error:
OGR: Unsupported geometry type
Error in CPL_sfc_from_wkt(x) : OGR error

Do I need to reformat the .geo variable first? If yes, how?


texts: QGIS 3.10 Issue on Export Selected Feature

texts:
I recently come across a weird issue when performing "Select by location" and then "Exporting selected features".

Basically, I have a large dataset (points) and a target area (polygon). I just need to select all points within the polygon. Sounds pretty easy and straight forward.

Firstly, I run "Select by location" and get around 74K features selected.
Secondly, I export those points by "right click -> export -> save selected features as".
Then, the exporting task will run for approx 15 mins but give me a complete empty layer!

I tried the same step on small datasets and the export works fine.
Does anyone have any advice?


texts:
I recently started using the GeoPackage format for vector data and I am very pleased with it. Its documentation mentions that it can also store raster data.

How can those be imported in an existing (or new) GeoPackage using QGIS or other open-source software, and what are its limitations?


texts:
I am using django-cities and geopy. I am trying to get city name, id and country name, id based on given latitude and longitude.

Normally this is how you can get city and country using geopy.

from geopy.geocoders import Nominatim

geolocator = Nominatim(timeout=3)
geolocator.reverse('52.5094982,13.3765983')
loc = location.raw
loc_dict = location.raw
print(loc_dict['address'])
and this is the output:

{'suburb': 'Tiergarten', 'country_code': 'de', 'country': 'Deutschland', 'postcode': '10117', 'attraction': 'Potsdamer Platz', 'road': 'Potsdamer Platz', 'city': 'Berlin', 'city_district': 'Mitte', 'state': 'Berlin'}

In django this is how i'm trying to get them:

Normally it should work, my problem is that geopy returns country: Deutschland whereas in django-cities it's Germany. And this happens to cities aswell, geopy gives me București and django-cities Bucharest.

How can i get city name and country name like in django-cities? If i get the name like in django-cities i can get the id.


texts: Leaflet check if L.geoJSON layer has any features

texts:
I am making my own layer control in a way.
I want to check if a named L.geoJSON has any features. If not I will fetch GeoJSON from an API and add it to the layer.

But I cannot work out how to check if it has any feature?
Have tried to find someway around it, have been trying to look for properties on the object, but cannot get it to work.

const mygeojsonLayer= L.geoJSON(null)

if(mygeojsonLayer._layers.length === 0){
        mygeojsonLayer.addTo(map)
} else {fetchData()}


texts: What does 'Weight Field' in IDW interpolation mean?

texts:
I am trying to figure out how the Weight Field option of Geostatistical Analyst IDW in ArcMap 10.3 works.

On IDW - ArcGIS Desktop Help I can find just this simple explanation:

  Used to emphasize an observation. The larger the weight, the more impact it has on the prediction. For coincident observations, assign the largest weight to the most reliable measurement.

What is the algorithm that creates the weight which emphasizes the observation?


texts: Loading vector layers from database/ schema using login username and password in QGIS?

texts:
I have a QGIS setup with a PostgreSQL / PostGIS database on a server and multiple clients.

I have created schemas in the database for each user. With a username password combination, each user is able to connect to his/ her schema to load his/ her vector layers.

I want to create a common QGIS project for all users which users will download onto the client PCs from a shared location. When the project file is opened, the users should enter their username and password to get connected to their respective schema, and the layers should get loaded.

I am able to do this by creating separate project files for each user, but how do I achieve it with a single project file?


texts:
I have a shapefile with, for example, X_coords and Y_coords of 870058.5631 and 3226158.25 respectively. I want to convert these into decimal degrees in ArcGIS pro, but none of the input coordinate formats seem to work.
They are along the Louisiana Coastline. The spatial reference is NAD 1983 UTM Zone 15N.


texts: Averaging set of rasters using Raster Calculator of ArcGIS Desktop?

texts:
How do I average a set of rasters?

I have remotely sensed chlorophyll a rasters for several days and I need to get an average of the days represented as one raster layer.

I am using ArcMap 9.2.


texts: Why do I get an error code 503 response when trying to access my CartoDB dashboard?

texts:
One pressing problem is that when I try to access my tables or visualizations on the CartoDB web app, I get this error page. My account is 'safedev'

Something that may or may not be relevant, but was not true before now:

When I try to q=SELECT CDB_UserTable()&api_key=[my-api-key] with the SQL api (via POST, but not via GET), I get:

          503 Service Unavailable           Error 503 Service Unavailable     Service Unavailable     Guru Meditation:     XID: 175897971 (EDIT: Will not match this number exactly - Jake)          Varnish cache server

However, if I do not supply my api_key (ie. q=SELECT CDB_UserTables()), I get normal results.

Also I can still get data from my tables via the SQL API, with or without api_key. For example,
q=SELECT * FROM "table"&api_key=[my-api-key] returns what I expect.


texts: Plot a discrete raster over a hillshade using ggplot2

texts:
I'm having a hard time fixing a somwhat simple ggplot. I have a raster of a hillshade and I would like to overlay it with a classified prediction raster with the viridis color scheme.
My code at the moment looks like this:
  # hillshade dataframe
  df.hsd = rasterToPoints(cropped_hillshades[[i]]) %>% as.data.frame(.)
  b.hs = seq(min(df.hsd[[3]]),max(df.hsd[[3]]),length.out=100)

  plot = ggplot(data=df.hsd, aes(x=x, y=y)) +
    geom_raster(aes(fill=df.hsd[[3]])) +
    scale_fill_gradientn(colours = grey(1:100/100), breaks=b.hs, guide="none") +
    geom_raster(data=df.pred, aes(x=x, y=y, fill=.data[["class"]]), alpha=.8) +
    scale_fill_viridis_c(name="Class") +
    coord_equal() +
    theme_minimal()

When I plot only the hillshade I get somehting like this:

When I only plot the susceptibility raster I get:

Both kind of look ok, but when I execute the lines above I get:
Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

And that's not really what I want as the Hillshade isn't visible anymore... I'm really done with trial and error in ggplot and would be grateful for any help.


texts: Converting multipe rasters to polygon through iteration fails

texts:
I need to convert several Rasterfiles (TIF) to polygon. To give some context: I saved the TIF-files in a local folder and I want to store the converted polygons in a database (to keep the information of shape length, etc.).

It works if I use the tool from modelbuilder: Raster to Polygon, for a single TIF file. But if I want to use the iterator; it always fails at the second iteration.

I make use of inline variables in the output of Raster to Polygon, but the problem is that this tool automatically adds an extension to the name (.shp) , which lead to a fail at the second iteration (Extension is manually removed at first iteration)

Does anyone has a solution for this?

I have to merge afterwards all this files to 1 polygonfile, I tried to other way around (Mosaic to new Raster and then Raster to Polygon ==> but this takes ages)

MODEL

GeoTiff = workspace (folder on C-drive) : all files start with Grnkrt_* / The output is saved to my default.gdb
FYI: I tried all possible inline variables already (%n%, %name,..)


texts:
I'm using GeoServer + PostGIS and I got an image generation error when I test my layer in WMS GetMap request:

http://localhost:8089/geoserver/pmc/wms?service=WMS&version=1.1.0&request=GetMap&layers=pmc:utbs&styles=&bbox=272206.156,7448903.5,303656.5,7484798.0&width=448&height=512&srs=EPSG:4326&format=image/jpeg

This is the image error:


texts: Mouse wheel zoom on leaflet map only after click

texts:
I want to make the zoom on mouse wheel only active after clicking once on the map, as described in Leaflet mouse wheel zoom only after click on map, but as a newby to all of this, I don't know where to put the code posted in the above mentioned thread.

This is my code so far:


texts: Is it possible to install Geoserver on Zentyal?

texts:
We want to setup a small WebGIS environment based on Geoserver/PostGIS or even the OpenGeoSuite if possible.
Our Server runs on Zentyal. Is it possible to install  Geoserver and PostGIS on Zentyal and if so is there any documentation available?
Or what would be an easy to operate platform so that we can create our WebGIS (we do not have much experience in running web servers though).


texts: Automate conversion from SHP to GPX file with Python

texts:
I've been trying to establish a method whereby I can automatically convert a shapefile to a GPX file using python, but having terrible luck in doing so.

Currently the only method that i've come across is using ogr2ogr (GDAL library) and then parsing that through os.system.

os.system(('ogr2ogr -f GPX %s %s') % ((os.path.join(toDirectory,nameOfGPX +".gpx"), inputLineShpfile)))

Though the end goal is to get a product of how the DNR Garmin spits out a GPX (see below).

As you see below it doesn't look good with multipart polygons looking like a mess.

ogr2ogr using GDAL:

End goal (using a manual process of DNR Garmin to convert)

Has anyone found a successful way to convert to gpx? I've got access to a basic version of ArcGIS and open to any open source or python libraries to assist.


texts: Access actions of ToolbarMenu/PanelMenu using PyQGIS

texts:
For a Python plugin I need to remove/disable some UI elements so that they can't be accessed by the user. Thus, I removed the respective actions, toolbars and menus from the interface.
However, in order to prevent the user from turning back on the disabled elements, their corresponding menu entries also have to be removed. Unfortunately, I didn't find a way to access the viewmenu's entries "Toolbars" and "Panels" (which contain the desired actions). They are not returned when calling the viewmenu's children() method. Likewise, the findChild()/findChildren() methods do not return anything helpful either.
I already had a look on the QGIS git repository and it seems that the desired menus and actions are not added to the viewmenu as its child elements. They seem to do not have any parent element at all..
So does anybody know how I can access the "Panels" and "Toolbars" menus from out my plugin code?


texts: Counting number of features in legend for each page when using Data Driven Pages?

texts: Is there a way to count number of features for each data driven and reflect it in the legend?


texts: Creating a DEM from irregularly spaced contour lines in QGIS

texts: I'm using QGIS to work with a contour line shapefile with contour lines that are not spaced evenly (e.g. not just 100 meter elevation increase between each line). The increase ranges from a 100 meter increase to a 5000 meter increase. Is there anyway to create a DEM from this data, or do I need to find a way to regularly space the elevation data?


texts: Why is my Image Mosaic in ArcGIS for Desktop showing a triangle?

texts:
When I mosaic two rasters in ArcGIS 10.2 for Desktop using Mosaic (Data Management) tool they give this image as shown below.

Kindly tell me what am I doing wrong?


texts: List the projection & datum of all the files in a folder

texts: In QGIS/GRASS, is it possible to list the projection details of all the shp,tab,dxf files in a folder ? We are having 100 - 150 files in too many folders. We would like to create a table of Proj details of all the files in those folders. Any idea ?


texts: Converting sqlite (spatialite) data to oracle sdo_geometry (using .NET)

texts:
Using .NET, I have a sqlite file that I have to import to oracle.  The geometry needs to be converted from sqlite (spatialite) geometry column to a sdo_geometry column in oracle:

// Loop through all records in sqlite file
foreach (DataRow row in dsSqLite.Tables["Default"].Rows)
{
  // Get geometry field from sqlite and get it ready for oracle format
  NetSdoGeometry.sdogeometry geom = row["GEOMETRY"];  // ERROR HERE
}

When I execute the above, I get an error "Unable to cast object of type 'System.Byte[]' to type 'NetSdoGeometry.sdogeometry'.  Also I do not see where to specify the SRID.  Does anyone know the correct way to convert the geometry column from spatialite to oracle sdo_geometry?

References:
NetSdoGeometry library: https://github.com/mapspiral/NetSdoGeometry


texts: how does one reset map.fitbounds (zoomtofeature()) to the original zoom?

texts:
I have a lovely map setup using Leaflet with MapBox tiles (only).  It is a United States map with polygons for each state with an onclick event to zoomToFeature.

The basic zoomToFeature() function from the tutorials does exactly what I want on the way in, but is there a similarly easy way to zoom out to original zoom?

I guess I could set a global variable to hold the state on the first event, and do an if inside zoomToFeature to either zoom or pull back, I just wonder if there is already some code that does this more elegantly.

Specifically for it to zoom out to exactly where the map was.  Not only to the original zoom level, but also the original center.


texts: Need assistance in project regarding Change Detection

texts:
My task initially was to classify the area into land-use classes like forest, sparse vegetation, agriculture, river, etc
I had done that by clustering the whole area and then reclassing it.
I had to further show some analysis and accuracy of my classification so that the change detection can be visualized.
So I need some guidance on how to approach next as I am stuck here and don’t know how to proceed further.

I am sending my code snippet of the year 2000.

I have data from 1995 to 2015 with 5 years interval

var poi=ee.Geometry.Point([78.032188,30.316496]).buffer(35000);
var dataset = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
              .filterDate('2000-01-01', '2000-12-31');
var trueColor321 = dataset.select(['B3', 'B2', 'B1']);
var tru=trueColor321.reduce(ee.Reducer.mean());
var truclip= tru.clip(poi);
var trueColor321Vis = {
       min: 0.0,
       max: 0.4,
       gamma: 1.2,
};
Map.setCenter(78.032188,30.316496, 9);
Map.addLayer(truclip, trueColor321Vis, 'True Color (321)');

var dataset1=ee.ImageCollection('LANDSAT/LT05/C01/T1_32DAY_NDVI').filterDate('2000-01-01','2000-12-31');
var NDVI= dataset1.select('NDVI');
var ndvimean=NDVI.reduce(ee.Reducer.mean());
var ndviclip=ndvimean.clip(poi);
var visparam={
   min: 0.0,
   max: 1.0,
   palette: [
'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
'66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
'012E01', '011D01', '011301'
     ]
    };
Map.addLayer(ndviclip,visparam,'NDVI',0);

//clustering
var training=ndviclip.sample({
region : poi,
scale : 30,
numPixels : 5000
 });
var clusterer=ee.Clusterer.wekaKMeans(20).train(training);
var result = ndviclip.cluster(clusterer);

Map.addLayer(result.randomVisualizer(),{},"clusters");

//forest classification
var reclass= result.remap([8,10,13,6,11,17,19],[1,1,1,1,1,1,1]);
Map.addLayer(reclass,{},'forest',0);

//water bodies
var waterbodies=result.remap([2],[2]);
Map.addLayer(waterbodies,{},'water',0);

//urban/barren land
var urban=result.remap([7,12,3,5,0,16,15],[3,3,3,3,3,3,3]);
Map.addLayer(urban,{},'urban/barren');

//agricultural
var agriculture=result.remap([4,8,17,19,18],[4,4,4,4,4]);
Map.addLayer(agriculture,{},'agriculture');

//sparse
var sparse= result.remap([1,9,4,14],[5,5,5,5]);
Map.addLayer(sparse,{},'sparse');


texts: How to create pre-configured tags for new ASP.NET Web Applications?

texts:
I'm using the "ArcGIS Server Manager" to create "Web Applications" with my maps. In ASP.NET project has a "web.config". How do I get my "web.config" come with some "tags" standard?

More specifically, these tags would be:

    < configuration> < configSections> < /configSections> < system.net> < settings> < >>servicePointManager expect100Continue="false" /> < /settings> < /system.net> < >>/configuration>

These tags are described in:
How do I use a layer of Bing / Google a "WebApplication" in "ArcGIS Server Manager"?


texts: What is the relationship between SDSFIE and C/JMTK?

texts:
C/JMTK = Commercial Joint Mapping Toolkit  = http://www.cjmtk.com/
SDSFIE = Spatial Data Standard for Facilities, Infrastructure & Environment = http://www.sdsfie.org

I've been working on a project using the C/JMTK, so I'm pretty familiar with it. It is basically a repacked version of ESRI's ArcGIS and ArcObjects technologies. Recently I was asked if our databases meet the SDSFIE standard, or could we work with SDSFIE-standard databases.

My question is what is the relationship between SDSFIE and C/JMTK? It seems to me after some preliminary reading that SDSFIE is a description of how the database schema is developed. Is that correct? Is it just schema or are there other requirements that restrict the underlying technology choices (e.g. geodatabase format or sql database vendor).

Can the C/JMTK (e.g. ESRI ArcGIS & ArcObjects technologies) support the SDSFIE standard?

Thanks!

P.S. Tried to tag this with C/JMTK and SDSFIE but I don't have enough reputation points.


texts: How to access script (source) of custom tool stored in tbx?

texts:
I have few custom toolboxes, with script tools.
I can access individual toolbar by:

ToolBar=arcpy.ImportToolbox(toolbar_path)

I can access each tool by iterating through

ToolBar.__all__

I’d like to search/find key word inside tool’s Python script.

How can I access it?


texts: How to automatically plot intermediate contours with QGIS, GRASS GIS or SAGA GIS?

texts:
I'm mapping a steep mountain range. My main contour interval is 50 m and my intermediate contour interval is 10 m. In steep areas, the interval of 10 m is too small, so I need to leave of the 10 m-interval and use 20 m or 25 m instead.

Is there a way, to do that automated in QGIS, SAGA GIS or GRASS GIS?


texts: How to merge USGS contour data from files into state?

texts:
Long story short, I'm looking to take 1/3 arc-second 1x1 degree contour shapefiles from the USGS National Map Viewer and transform that into a contoured map of a particular state, but am encountering some significant processing problems along the way. Question is in bold near the end of the post.

To download the maps, I selected all 3DEP Elevation Products, then downloaded the 1/3 arc-second 1x1 degree quads that covered at least part of the state I'm interested in (California). This was easy.

I then used ogr2ogr to merge the contour quads into a large grid with the help of a shell script that I wrote, shown below:

#!/bin/bash

ogr2ogrdir=''
directory='/'
output=''

cd $directory || exit

i=1

while IFS= read -r -d '' shpfile
do
    if [ ! -e $directory/$output ]
        then
            echo "[$i] Creating $output from shapefile $shpfile..."
            $ogr2ogrdir -f 'esri shapefile' $directory/$output "$shpfile"
        else
            echo "[$i] Appending shapefile $shpfile to $output..."
            $ogr2ogrdir -f 'esri shapefile' -update -append $directory/$output "$shpfile"
    fi
    ((i++))
done < <(find . -name '*.shp' -print0)

However, naturally ogr2ogr does not support the generation of files larger than 4GB. California was looking likely to be closer to 6GB. To work around this, I picked an arbitrary minute of latitude to divide California into northern and southern sections, and then merged the maps into two files. This produced two shapefiles of approximately 3GB each, looking like the following images. After some processing to remove contour lines outside of the state boundary (Vector > Geoprocessing Tools > Clip using a state boundary polygon as the clip mask), here is what the two files look like from a screenshot from QGIS below.

What I'd now like to do is merge these two separate files into a single state shapefile of contour lines. But, I'm encountering issues with QGIS. I've tried Vector > Data Management Tools > Merge Vector Layers in QGIS, but sadly the output is not what I'm expecting. As you can see in the image below, some (why?) of the southern CA original quads in the southern CA shapefile are missing in the final merged product.

I'm not even sure why the artifacting is even occurring along the original quadrangle lines - ogr2ogr supposedly merged them, or at least that's my understanding of the situation.

How can I resolve this issue to produce a clean output of contour lines for the state of California in a single shapefile?

Additionally, I'm noticing at the boundary between my northern and southern demarcation point, I'm getting some misalignment between the contour files. Shown in the image below is the contour lines intersecting at the western point of the state at 41 degrees latitude. Could this be causing the error I'm seeing above?


texts: Generate points inside multipolygons using GeoPandas

texts:
I have a shapefile which has multipolygon geometries. I want to generate random points inside them using GeoPandas. I tried with the below code but it's not giving geodataframe of points.
df = gpd.read_file(multipolygonshp)
geom= df['geometry']
final_geom = geom.unary_union
x_min, y_min, x_max, y_max = final_geom.bounds

n = 3000
x = np.random.uniform(x_min, x_max, n)
y = np.random.uniform(y_min, y_max, n)

gdf_points = gpd.GeoSeries(gpd.points_from_xy(x, y))
gdf_points = gpd.GeoSeries(gdf_points.within(final_geom.unary_union))

<ipython-input-125-ed976bd8ef24>:23: FutureWarning:     You are passing non-geometry data to the GeoSeries constructor. Currently,
    it falls back to returning a pandas Series. But in the future, we will start
    to raise a TypeError instead.
  gdf_points = gpd.GeoSeries(gdf_points.within(final_geom))


texts: Publishing ArcGIS Desktop 10.0 Custom Toolbox to ArcGIS Server 10.6?

texts:
I am pretty new to this stuff.

I have a Custom Toolbox with python scripts that was created in ArcGIS Desktop 10.0.

I am upgrading to ArcGIS Desktop and Server 10.6 and trying to publish the Toolbox as a service.

I need a .sd file created but seem to hit a dead end on finding out how to do this.

I am able to open the .tbx in ArcCatalog but nothing shows up in the Results.

What am I missing?


texts:
user borrowed from friend

since i use vector layer with wfs protocol, there is a little while before features loaded on the map. So i need to tell user to wait a moment until they load.

what best way to do that ?
i found this http://osgeo-org.1560.n6.nabble.com/Patch-quot-waiting-quot-mouse-cursor-on-WFS-loading-td3961983.html

but not sure if it will work. any body tried it ?


texts: Creating n randomly distributed points within given (multi)polygon with PyQGIS?

texts:
Is there a way to create n randomly distributed point features within a given polygon using QGIS, PyQGIS and/or PostGIS?

My first approach was to implement something like (fantasy language)

n = 0
repeat until n = 20
    create point feature within polygons boundingbox
    if point within polgon
        increment n
        keep point

It's easy to create randomly distributed points within a rectangle, but depending on the shape of the polygon, this algorithm may create a huge amount of points that are not within the polygon and may run a while, if it terminates at all...


texts: Accessing arcpy.mapping .tif layer object?

texts:
Can someone tell me why I am getting this value error when attempting to access the Montreal.tif layer object?

The first three lines access the object and set the transparency perfectly well, but when I attempt to access it using the following 4 lines I get the illustrated ValueError?

I need the second method since the b variable (Montreal) will change.

ArcGIS Desktop v10.1, ArcInfo license.

The following is a successful workaround (among others), but I'm still not sure why the above was giving me an error?

mxd = arcpy.mapping.MapDocument("Current")
df = arcpy.mapping.ListDataFrames(mxd)[0]

b = "Montreal"
for lyr in arcpy.mapping.ListLayers(mxd, "", df):
    if "{0}.tif".format(b) in lyr.name:
        lyr.transparency = 50


texts: Convex Hull not outputting expected number of features

texts:
I have a .csv file with observations of various species and I want a polygon to illustrate the area where those species occur. I ran the Convex Hull, grouped by the species' names field, which returned a feature for each species, as expected.

The problem is that there were some species that weren't expected to return an error because they don't have enough geometries to form a polygon, but the number of errors was higher than the number of species with 1 or 2 observations.

What causes the Convex Hull in QGIS 3.14 to return an error?


texts: Adding fields in ArcMap to enable entering of range numbers?

texts: I'm trying to enter range numbers such as 1 - 5 , 4- 40, etc in an attribute table for one of my shapefiles as I'm trying to display the depth ranges of different areas around my study area with different colours. But I'm not sure what 'Type' to choose from when adding fields or if this is even the right way to display this type of data. I'm only trying to input the numbers, not the words such as 'Shallow' etc that come with it.


texts:
I've digitised some symbols in QGIS (2.14.1) as a new Shapefile point layer. I want the symbols to be rotated to different angles.

I've worked out that I can have a data defined rotation field in my attribute table and I've put this in and it seems to be working. However I can't workout if there is a way to rotate the symbol in the map window and have that update in the attribute table: so I'm stuck with putting values into my attribute table manually.

I know how I can rotate label with the Label Toolbar so I really hope that there is a similar function for symbol rotation.


texts: Gettting all countries information in world using Google API

texts:
The Scenario:
My customer wants to let the users insert their location either by choosing the country, city, region, state, etc.. from dropdowns OR the user can choose the location using google API and then fill the dropdowns by the data comes from google API after comparing the retrieved values with the values in the dropdowns.
to accomplish this scenario I should save all countries Info in a local database to fill dropdowns and of course should be the same names in google API to match it later, and then the user will be able to select their location using dropdowns.
So, My question is it possible to get all countries Info from google API or any library based on google API to use it in the dropdowns.
I found this answer before but the answer is old and the source of data in Geonames not based on google API, So if I use it I will face a lot of differences in names of the cities, regions, states, etc..
and what I want is to get all countries data based on google API even if this service paid.
Additional Note: I Use ( C#, Asp.Net Core, Xamarin Forms, Angular )


texts:
I have a PointZM shapefile and I would like to drop the Z value and just write out the PointM values to a new shapefile (I'm trying to test a shapefile reader so other formats are not an option).

If I specify the geometry in the export window and uncheck include z-dimension I get a simple POINT file with just X & Y. If I check that box I get a POINTZ file with X, Y & Z.

Is there some trick I can use or should I file an enhancement request?


texts: Get average values in some specific range based on spatial analysis with QGIS

texts:
I'm working on QGIS to compute average values attached to polygons around a line.

What I'd like to do is compute the average values of polygon datawithin a user defined range from a line.
How would I go about doing this?

I have attached a picture below for reference:


texts: Python - Calculating Burned Area Index results in empty plot

texts:
I've followed this tutorial to calculate the dNBR Index, using the given sample tif files which contain a burned area. My results were as shown on the tutorial. So far so good.
Then I tried to calculate the Burned Area Index (BAI) after making a few modifications on the above tutorial's code. The result is an empty plot.
This is my implementation:
(Note: The code downloads tif images from here)
from glob import glob
import matplotlib.pyplot as plt
import xarray as xr
import rioxarray as rxr
import earthpy.plot as ep
import earthpy as et
import os

def combine_tifs(tif_list):
    """A function that combines a list of tifs in the same CRS
    and of the same extent into an xarray object

    Parameters
    ----------
    tif_list : list
    A list of paths to the tif files that you wish to combine.

    Returns
    -------
    An xarray object with all of the tif files in the listmerged into
    a single object.

    """

    out_xr = []
    for i, tif_path in enumerate(tif_list):
    out_xr.append(rxr.open_rasterio(tif_path, masked=True).squeeze())
    out_xr[i]["band"] = i+1

    return xr.concat(out_xr, dim="band")

path_to_downloaded_data = et.data.get_data('cold-springs-fire', replace=True)

all_landsat_bands_path = glob(os.path.join(path_to_downloaded_data,
                                       "landsat_collect",
                                       "LC080340322016072301T1-SC20180214145802",
                                       "crop",
                                       "*band[1-7]*.tif"))

all_landsat_bands_path.sort()
all_bands = combine_tifs(all_landsat_bands_path)

red_square = (0.1 - all_bands[3]) * (0.1 - all_bands[3])
nir_square = (0.06 - all_bands[4]) * (0.06 - all_bands[4])
burned_area_index = 1 / (red_square + nir_square)

fig1, ax = plt.subplots(figsize=(12, 12))
fig1 = plt.gcf()
ep.plot_bands(burned_area_index,
          cmap='PiYG',
          scale=False,
          cbar=False,  # Change this to TRUE to display the color bar
          ax=ax,
          vmin=-1, vmax=1,
          title="BAI")
plt.show()

The values of the burned_area_index array are very low, in the magnitude of 10^-7.
My question is, if the formula is correct, what might be the problem that results in an empty plot?


texts:
I've got a vector of points, polygons, and lines. And I'm trying to convert the points into a GML file like a Web Feature Service. I've been doing some reading but can't find any methods that can do this?

I'm coding in C++ with QT

Using GDAL Libraries


texts: Getting range value from exponential variogram model

texts:
I am using vgm() from the gstat library in R:

variogr <-vgm(sill,"Exp",range=range,nugget=0)
output <-krige(...,model=variogr)

As you can see I'm using fixed psill and range values.

Before this, another script (on which I have not control) is evaluating the variogram parameters. Namely, I know that my points fit well in a p1*(1-exp(-x/p2)), and I am given the exact values of p1 and p2. My question is how to get sill and range (sill and range) from the exponential model parameters p1 and p2?

I know that the sill value is == p1, because we are choosing nugget to be zero. But I'm having troubles finding which value should I give to the range. I have been told by my professor that range == p2, but from what I have read elsewhere there should also be a factor (like, p2 == k*range) that might have different values depending on the proper definition of range.

So, given p1 and p2, which value should I use as a range input to the vgm() function? Is it true that range == p2? What kind of range is vgm() expecting me to put in?


texts:
Is it a way to digitize indicating a curve (colored on a white or black or different colour background) from point to point ? I think i've been told that, but i don't KNOW to do it !
Please consider a quite low GIS knowleadge....(sorry)


texts:
I usually use the proportional circles extension of QGIS but I need to have proportional pork or pigs according to the number of pigs per exploitations (it's a map about agriculture).
How do I do that ?


texts: Labels still showing when no Label Features turned on in ArcMap?

texts: I'm continuing a project from someone else and even after I turn off all the layers, there are still labels displayed in the data frame. None of them have "label features" on, so I'm not sure where these labels are coming from. I'm using ArcMap 10.5.


texts:
I have a style pulled from my GeoServer and have created a basic legend, with a fixed PNG. However I was wondering if it is possible to have this legend be pulled from the layers SLD? so that if the SLD is changed the legend will update too?

This is the Legend code

    var legend = L.control({position: 'bottomright'});

    legend.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'info legend');

        div.innerHTML +=
        '<img src="SymbolDelete.png" alt="legend" width="134" height="147">';

    return div;
    };

    legend.addTo(map);


texts: Distance matrix tool, how to include distances at 0 metres?


texts: How to perform Random Forest land cover classification?

texts:
This is a follow-up to a previous post: Machine Learning Algorithms for Land Cover Classification.

It seems that the Random Forest (RF) classification method is gaining much momentum in the remote sensing world.  I am particularly interested in RF due to many of its strengths:

A nonparametric approach suited to remote sensing data
High reported classification accuracy
Variable importance is reported

Given these strengths, I would like to perform Random Forest land classification using high resolution 4 band imagery.  There is a lot of material and research touting the advantages of Random Forest, yet very little information exists on how to actually perform the classification analysis.  I am familiar with RF regression using R and would prefer to use this environment to run the RF classification algorithm.

How do I collect, process and input training data (i.e. based on high resolution CIR aerial imagery) into the Random Forest algorithm using R?  Any step-wise advice on how to produce a classified land cover raster would be greatly appreciated.


texts: How to check if a geoserver layer is raster or vector using gsconfig

texts:
I am using gsconfig and python in order to retrieve information regarding some geoserver layers.

What I am trying to do is to find out if a specific layer is raster or vector.
I read online that I need to use get_resource function to do that.

This is what I do:

name = "layer_name"
workspace = "my_workspace"
cat = Catalog(location, user,pass)
resource = cat.get_resource(name, workspace=workspace)

The question is how to get the type of the layer after that. When I print the resource I get:

 <geoserver.resource.Coverage object at 0x7fb2f8ffca90>

And if the layer is vector then I get:

 <geoserver.resource.FeatureType object at 0x7fb2f9e64c50>

What's the way to access the properties of the layer and see if its vector or raster so later I can do something like this:

if layer = "raster":
   # do something
elif layer = "vector":
  # do something else

Doing:

 print dir(resource)

I get back a list with the attributes of the class:

['__class__', '__delattr__', '__dict__', '__doc__', '__format__', '__getattribute__', '__hash__', '__init__', '__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_href', '_store', '_workspace', 'abstract', 'advertised', 'catalog', 'clear', 'dirty', 'dom', 'enabled', 'fetch', 'href', 'keywords', 'latlon_bbox', 'message', 'metadata', 'metadata_links', 'name', 'native_bbox', 'projection', 'projection_policy', 'refresh', 'request_srs_list', 'resource_type', 'response_srs_list', 'save_method', 'serialize', 'store', 'supported_formats', 'title', 'url_part_stores', 'url_part_types', 'workspace', 'writers']

With help from @iant, I found a way to see if my layer is vector or raster. I am not sure if this is the best way to do this but it works.

I actually followed step by step @iant's answer but in my if statements instead of putting:

if type(resource) is geoserver.resource.Coverage:
    # this wouldn't work properly.

I put this:

a = str(type(resource))
if a == "<class 'geoserver.resource.Coverage'>":

So I actually transformed the result of type(resource) into a string.
It should be the case that the first way does work. I am not sure why it fails.


texts: vuejs and openlayers, constructor is not a constructor

texts:
Uncaught TypeError: Constructor is not a constructor
    at Draw.geometryFunction [as geometryFunction_] (VM2096 app.js:40101)

Any idea what causes that? after that, each mouse movement makes an error.

I am trying to draw a polygon
(Laravel, VueJS, OpenLayers)


texts:
should my "site manager" url contain "localhost:6080"? or should it be "(my pc name):6080"? After giving arcgis server account permissions to proper folders I created my site, but it has the name "localhost:6080", will I still be able to link published services with REST Url when making web maps? I am on ags server 10.1.

I have a .com registered and I am making a subdomain at the moment. How do I point my arcgis server to this domain? (the process to making my files accessible outside my network) or do I need a webserver and arcgis for server installed to do what I wish to do?(making web maps accessible to people on the internet). I do not have a web server installed. I am trying to install arcgis for server 10.1 on my windows 7 ultimate machine.

According to Branden, using "http://gis.mycompanywebsite.com/arcgis/rest/services/101/my_service_name/MapServer" Will point my javascript to my services, using this will it matter if I have localhost in my url, if so, how do I avoid having a localhost site?


texts: Export to drive only exports black and white images

texts:
I am generating a burn severity map for the assessment of the areas affected by wildfires in California and Australia.using this code:
https://code.earthengine.google.com/b455ba8cf4b5bee822bb7ff8935e6209
My problem is while trying to export the image shown in earth engine editor. The Output is only in black and white.
Can someone help me ?
/==========================================================================================
//                                  PREPARE FILE EXPORT

var id = dNBR.id().getInfo();

Export.image.toDrive({image: dNBR, scale: 30, description: id, fileNamePrefix: 'dNBR',
  region: area, maxPixels: 1e10});


texts: Field mapping features to different feature with different field names using ArcPy

texts:
I was trying to you this instead of field mapping because field mapping is hard to understand.

outfc = r'\\water.com\Blank.gdb\SewerCollection\ssLateralLine'
infc = r'\\water.com\Test.gdb\PW_Waste_Water\WastewaterUtilitySegment'

#The position of the fields must match, Name=Description etc.
outfcfield = ['INSTALLDATE', 'FACILITYID', 'SIZE', 'InstallationID', 'DOWNELEV', 'UPELEV', 'SIZE', 'LIFECYCLE', 'MATERIAL', 'PROJECTID_ADD', 'SEGMENTTYPE']
infcfields = ['equipDatePutInService', 'featureName', 'insideDiameter', 'installationId', 'invertElevationDownstream', 'invertElevationUom', 'nominalPipeSize', 'operationalStatus', 'pipeMaterial', 'projectId', 'wastewaterSegmentType']
icur = arcpy.da.InsertCursor(infc, infcfields)
with arcpy.da.SearchCursor(outfc, outfcfields) as cursor:
        for row in cursor:
            try:
                icur.insertRow(row)
            except:
                pass
            print "row[0:]: "+str(row[0:])
del icur

The insert cursor script above work but it didn't add the features just there data. I need to copy my feature to a different GDB that has a different Schema.
I found the script below here https://community.esri.com/thread/185431-append-tool-and-field-mapping-help-examples and thought i could ad my data and see what happens then go from there.

import arcpy

append_layer = r'\\water.com\Test.gdb\PW_Waste_Water\WastewaterUtilitySegment'
target_layer = r'\\water.com\Blank.gdb\SewerCollection\ssLateralLine'

#This object looks like the empty grid of fields
#you see when you first open the append tool in the toolbox
fieldmappings = arcpy.FieldMappings()

#Like when you manually choose a layer in the toolbox and it adds the fields to grid
fieldmappings.addTable(target_layer)
fieldmappings.addTable(append_layer)

#####Lets map fields that have different names!
list_of_fields_we_will_map = []
#Lets chuck some tuples into the list we made
list_of_fields_we_will_map.append(('equipDatePutInService', 'facilityNumber', 'installationId', 'numberOfPumps', 'operationalStatus',   'projectId','featureName'))

for field_map in list_of_fields_we_will_map:
    #Find the fields index by name. e.g 'TaxPin'
    field_to_map_index = fieldmappings.findFieldMapIndex(field_map[0])
    #Grab "A copy" of the current field map object for this particular field
    field_to_map = fieldmappings.getFieldMap(field_to_map_index)
    #Update its data source to add the input from the the append layer
    field_to_map.addInputField(append_layer, field_map[1])
    #We edited a copy, update our data grid object with it
    fieldmappings.replaceFieldMap(field_to_map_index, field_to_map)

#Create a list of append datasets and run the the tool
inData = [append_layer]
arcpy.Append_management(inData, target_layer, field_mappings=fieldmappings)

Runtime error
Traceback (most recent call last):
  File "<string>", line 25, in <module>
  File "d:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\arcobjects\arcobjects.py", line 328, in getFieldMap
return convertArcObjectToPythonObject(self._arc_object.GetFieldMap(*gp_fixargs(args)))

RuntimeError: FieldMappings: Error in getting field map from field mapping for GetFieldMap

How do I field map my ssLateralLine feature to the WastewaterUtilitySegment features?


texts: Amnesic installation of QGIS 3.4 on Ubuntu 18.04

texts:
I installed QGIS 3.4 with the repository https://qgis.org/ubuntu-ltr bionic main on an Ubuntu 18.04 machine.
I had some troubles due to the presence of several repositories listed in my file sources.list. But now the software is installed and is properly running.

However its settings are not stable : neither the toolbar's display nor the location of the panels
nor the import of my personalized color palette nor even the recently opened project files
nor the connection to the databases are kept by the software.

I checked if the profile folder existed : /home/user/.local/share/QGIS/QGIS3/profiles/default.
It's here, no problem. It has even the correct writing accesses for the owner and its group (and is "read-only" for the other users).

Do you have some clues?


texts: GDAL scripts not working on QGIS3 High Sierra

texts:
Good day everyone,

I had the same issue that has already been posted with running GDAL scripts on MacOs with QGIS3. I am trying to clip a raster and used to get this error:

GDAL command output: /bin/sh: (scriptname).py: command not found

Then I followed the workaround from the readme.txt to use customs variable and now I have another error...

Traceback (most recent call last): File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/gdal/GdalAlgorithm.py", line 115, in processAlgorithm GdalUtils.runGdal(commands, feedback)
File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/gdal/GdalUtils.py", line 83, in runGdal if not path.lower() in envval.lower().split(os.pathsep):
AttributeError: 'NoneType' object has no attribute 'lower'

I also tried to re-install GDAL after QGIS3 but no luck...

I am running High Sierra 10.13.4 and QGIS3 3.0.2.
For info, it doesn't work on QGIS 2 either...

What can I do?


texts: Converting LAS data from buildings to detailed shapes using ArcGIS Pro?

texts:
What we want to achieve is to convert LAS data to a GIS format like shapes. On the internet, there are enough examples how to convert the LAS data height to a building footprint.

What we want is extra detail, the basic shape of the footprint should be split where the building changes in height.

How do we convert LAS data to a footprint shape where (with the right settings) it splits a shape into more shapes with height like this?

Example of what to achieve:

It would be nice if the multipatch feature would allow this, but so far we didn't find an answer in that direction.


texts: GeoDjango - Find all locations within X miles of provided point using bounding box

texts:
I'm trying to implement some search functionality for a web app, and have run into an issue where I don't know the proper functions to execute this.

Given a Point, a bounding_box database field with the EWKB representation of the geometry and a given distance value (aka 10km), I want to find all records that are within the radius or intersect with the outer boundary.

An example bounding box is:

0103000020E61000000100000005000000800F04E8AE1751C09EEE3F0181234640F04A8F5A600E51C09EEE3F0181234640F04A8F5A600E51C079B75A336B314640800F04E8AE1751C079B75A336B314640800F04E8AE1751C09EEE3F0181234640

Given that this function will be executed many times, with possible caching, I'd like to avoid any ST_BUFFER usage.

There are multiple things happening here:

How can I create a bounding box representation using the EWKB value in the bounding_box field?  Aggregate functions like ST_POLOGONIZE are not allowed within WHERE blocks.
Is ST_DWITHIN usage appropriate here?
This function will be executed many times.  What sort of optimizations are suggested?


texts: Arcmap refuses to change my layout view to landscape when I change the options in Page and Print setup. It keeps the exact same portrait box as before.


texts: Why does the commands listgeo and gdalinfo don't show the same Proj4 Strings?

texts:
As the title says, why doesn't the unix tools listgeo and gdalinfo show the same proj4 strings?

Check this file:
https://github.com/geotrellis/geotrellis/blob/master/raster-test/data/econic.tif?raw=true

This gives with listgeo -proj4 econic.tif the proj4 string:

+proj=eqdc +lat_1=33.903634028 +lat_2=33.625290028 +lat_0=33.764462028 +lon_0=-117.474542889 +x_0=0.000 +y_0=0.000 +ellps=clrk66 +units=m

And with gdalinfo -proj4 econic.tif the proj4 string:

+proj=eqdc +lat_0=33.76446202777777 +lon_0=-117.4745428888889 +lat_1=33.90363402777778 +lat_2=33.62529002777778 +x_0=0 +y_0=0 +datum=NAD27 +units=m +no_defs

The only thing different is the ellps flag for listgeo and datum for gdalinfo. I have looked at the listgeo source and it never writes a datum flag.

Another question here of course is how do I compare these two, in a module I'm writing, and still get that they are the same?


texts: How to compare a set of images with a shapefile?

texts:
I have around 150 images (rasters), each one is just a map of certain "classes", e.g. pixel value 1=forest.

I also have a shapefile that consists of roughly the same area as the 150 images constitute when put together. The shapefile also includes the clasee 1=forest. How could I go around to compare if the "classes" match, and how well do they.

Any ideas?

I'm using Python and QGIS.


texts: Changing display name of layer with PyQGIS?

texts:
I did a batch processing of several vector layers and now the display name for every file is "dissolved" when accessing the layer.name() function in QGIS. I read in this answer from underdark to Changing layer name of output vector from processing script in QGIS?, that I can set the name to the file name but unfortunaly too late.

I am able to get all the files in my project to be in a list with

layers = iface.mapCanvas().layers()

and all the names with this:

names = [layer.name() for layer in QgsMapLayerRegistry.instance().mapLayers().values()]

With the help of this answer to Getting path of project, or layer file in PyQGIS?, I can get the names of the file I am working with:

import os
real_names = []
for li in layers:
   (myDirectory,nameFile) = os.path.split(li.dataProvider().dataSourceUri())
   real_names.append(nameFile.split("|")[0]

Reading the API docu for the QgsVectorLayer didnt helped me.
What I  want  is to rename the files to show the filename as display name. How to archive that with the python console?


texts:
As my project requirement I want to draw a polygon look like this (>. In this case draw polygon sides too dense for arc. Successfully draw this using LineString but how to draw it using polygon. Need polygon because I have to fill color in between.

code in OpenLayers 2

var linearRing = new OpenLayers.Geometry.LinearRing(vertices);
        return new OpenLayers.Geometry.Polygon([linearRing]);

Here vertices is the array of points successfully drawn like the image
In OpenLayers 3

var layerLines = new ol.layer.Vector({
        source: new ol.source.Vector({
            features: [new ol.Feature({
                geometry: new ol.geom.LineString(markers, 'XY'),

            })]
        })
        //style: iconStyle
    });

    map.addLayer(layerLines);

Working but can't fill the in between area.

If we replace geometry: new ol.geom.LineString(markers, 'XY'), line by this

geometry: new ol.geom.Polygon(markers, 'XY'),

then it can't draw a polygon


texts: How to join polygons and lines vector maps from openstreetmap?

texts:
I am struggling with the data format OSM gives me: I want to use OSM data for routing, but all round-abouts and parking lots are not in the lines layer. They are situated in the poly(lines) layer.

Thus I could not find a working way to combine both layers by their same type (e.g. highway=secondary) to have the roundabouts in place for valid routing.

I imagine that this is quite a easy one, but I am stuck. can someone help me out?

Frank.

After v.patch I added my stations to the network and added missing nodes

v.net --overwrite --verbose input=network points=stations output=network_tmp operation=connect thresh=500
v.net input=network_tmp operation=nodes output=network

BUT if I now want to do routing, I still get a lot of errors saying that the points are not reachable from each other. On the network itself it looks ok.

echo 2 2 4 |v.net.path --overwrite --verbose input=network output=path_2_4

I tried a cleanup as found on the web

v.clean input=network output=network_cleaned tool=break,rmdupl

This does not improve the routing. As the manual says it breaks the lines at each line intersection and therefore all information regarding bridges, tunnels etc. are lost.

Any idea why it is that hard to use openstreetmap data or are my steps wrong?


texts: Calculating mean slope: Harmonic or Arithmetic Mean?

texts:
I have to calculate the average upslope percent slope for a large dataset, the basic method is detailed here.   However, I've begun to wonder if harmonic mean might be more appropriate than the standard arithmetic mean, since it is technically a rate of change.  I haven't seen this turn up in any of the other discussions on averaging slope over points, areas, lines, etc.  It should be fairly straightforward to accomplish.

edit: The purpose of calculating average slope in this case is to generate one parameter (of many) to be used in modeling channel initiation thresholds.  I have a set of field-collected channel head locations that I will collect the flow accumulation, various average upslope parameters, etc, and will use multiple linear regression to try to describe accumulation thresholds in terms of the other parameters.


texts: PostGIS query for selecting features not connected to rest of road network?

texts:
I have a dataset of roads in a city stored in a PostGIS enabled database.  I would like to identify and remove any feature which is not connected to the rest of the road network.  This behavior is common around the edges of networks, where the data was clipped.

I have found other answers on the GIS SE that do this, but they all appear to use ArcGIS, and I would like to solve the problem using PostGIS, QGIS, or a Python script.

My datasets can include over a million features, so performance is a concern.


texts: How to attach pictures saved in relative path

texts:
I have a QGIS project saved under G:\USER\xyz\projects\QGIS.qgz

In the same folder is a subfolder Pictures. My layer contains a column with the picture names (For example: 101.jpg). The picture 101.jpg is within this subfolder Pictures. How can I tell QGIS to connect this picture name in the layer's column to that actual picture in the subfolder?

In the screenshot you can see my configuration - I think it should work like this? Relative path checked and in the default path field I wrote the name of the subfolder. It is not working though. Any thoughts? I want to achieve that the whole project folder is portable so I can use it with QField.


texts: Polygonize every single pixel of a raster

texts:
How can I get every pixel separated? When I calc unique values, if two pixels have the same value, it converts in one row in the final shape. I need every single pixel converted into one unique geometry.

Any ideas?

My code is:

def geojsonizer(ndvi_path, downloads_path, ndvi_date, UHC_name):
entry_date = datetime.today().strftime('%Y%m%d')
# Read input band with Rasterio
with rasterio.open(ndvi_path) as raster_ndvi:
    crs = raster_ndvi.crs
    huso = str(crs).replace(':', '')
    src_band = raster_ndvi.read(1)
    # Keep track of unique pixel values in the input band
    unique_values = np.unique(src_band)
    # Polygonize with Rasterio. `shapes()` returns an iterable
    # of (geom, value) as tuples
    shapes = list(rasterio.features.shapes(src_band, transform=raster_ndvi.transform))

geojson_schema = {
    'geometry': 'MultiPolygon',
    'properties': {'name': 'str', 'pixel_id': 'int', 'ndvi_value': 'float', 'ndvi_date': 'str',
                   'entry_date': 'str'}
}
with fiona.open(downloads_path + '\\' + UHC_name + '_' + huso + '_ndvi.geojson', 'w', driver='GeoJSON',
                schema=geojson_schema, crs=crs) as geojson:
    for pixel_value in unique_values:
        if pixel_value != 999 and pixel_value != 0:
            polygons = [shape(geom) for geom, value in shapes
                        if value == pixel_value]
            multipolygon = MultiPolygon(polygons)
            pixel_centroid = (str(multipolygon.centroid).strip('POINT ()')).replace(" ", "")
            geojson.write({
                'geometry': mapping(multipolygon),
                'properties': {'name': 'nombre_finca',
                               'pixel_id': pixel_centroid,
                               'ndvi_value': float(pixel_value),
                               'ndvi_date': str(ndvi_date),
                               'entry_date': entry_date,
                               }


texts: Get X, Y coordinates in the center of a raster cell

texts:
I am trying to get the X, Y coordinates of all the cells of a raster file, using Python 2.6 and ArcGIS 10.0.

I would like to get the X, Y coordinates in the center of the cell.

Is there any way to do so?


texts: I develop an application and I want to erase overlay features in a feature class in the application. I use ArcEngine 10. I don't have another layer that have erase features then I can't use erase tool in overlay toolset. How to erase using feature(s) within the same feature class? is there any tools or programmatically way to do this?


texts: earthEngine expression and ternary operator: why does evaluate false value when always true?

texts:
Why does EE `.expression('A?B:C') take into account C if A is always true?

I am using conditional operations with the image.expression('condition ? true_action : false_action) as documented here. My understanding of the A ? B : C constructs is that it is a ternary operator, that says: if A is true, return B. Otherwise, if A is false, return C.

I expect hence that if is A is always true, C will never be evaluated. So if C returns a non-sensical answer (like log(-1)) that should not affect my result. This does not seem to be the case in earthEngine: if C is non-sensical, the overall result will be show up with Masked pixel. Indeed compare:

var expr_ok = ee.Image(3).expression(" 3 <4   ? 1 : 2"): returns always 1
var expr_absurd = ee.Image(3).expression(" 3 <4   ? 1 : asin(2)"): does not return anything (Masked)

Is my understanding of ternary operator wrong, or is that a bug/unexpected feature of EE? I am aware I could use the image.where() approach, but my question is specifically about the .expression(A?B:C) approach.

Minimal code

Here is a very simple case showing the issue (my real case is more complicated, where I want to evaluate asin(x) only if x in [-1, 1])).

// different approaches
var expr_ok = ee.Image(3).expression(" 3 <4   ? 1 : 2")
var expr_absurd = ee.Image(3).expression(" 3 <4   ? 1 : asin(2)")

var where_ok = ee.Image(3).where(ee.Image(3).lte(4), 1).where(ee.Image(3).gte(4), 2)
var where_absurd = ee.Image(3).where(ee.Image(3).lte(4), 1).where(ee.Image(3).gte(4), ee.Image(3).asin())

//visu
var pal = {min: 1, max: 2, palette: ['0000FF', '00FF00', 'FF0000']}
Map.addLayer(expr_ok, pal, 'expr_ok', false);
Map.addLayer(expr_absurd, pal, 'expr_absurd');

Map.addLayer(where_ok, pal, 'where_ok', false);
Map.addLayer(where_absurd, pal, 'where_absurd', false);


texts: How to do selection with a WFS source in OpenLayers?

texts: One of my students wants to style a WFS layer and have selections too. I got as far as this code http://ian01.geog.psu.edu/geoserver/www/wfs/wfs-sld2.html which styles the WFS (though with out the text labels) and while the selection control seems to be working (i.e. the box is blue not red) nothing is actually selected.


texts: Handle a great amount of Markers with OpenLayers

texts:
My general Problem is to maintain a lot of markers. For the beginning I want to display all Bus Stops of Berlin, Germany and surroundings as Markers. I already got that I have to use a strategy like Startegy.BBOX for performance issuses (my dataset contains nearly 13000 Bus Stops). I found an Example where somebody gets his marker positions from a databases, parses it to a giant string and uses it with Format.Text. I guess this would work for me too, but I still wonder if thats the best way to do it. What if my dataset would be much bigger (maybe to display all bus stops worldwide)? I guess putting the whole content of the tatabase to a giant string and use Format.Text woulden't be a good idea than.

What would be the best way to do it? Is there a simple way to only load marker positions when they are actually needed instead of preloading them all?


texts: How to fix a IndexError: string index out of range

texts:
I'm still working on my script for vineyards. I need to count row and vine in row.
So I code this :

#Nombre de rang
nbr = QgsExpression( "ID LIKE '%-1'" )
itnbr = layer.getFeatures( QgsFeatureRequest( nbr ) )
litnbr = list(itnbr)

print 'nombre de rang : ' + str(len(litnbr))

#Nombre de pied dans le rang
for i in range(len(litnbr)):
    rang = str(i+1)
    print 'rang : '+ rang
    itp = layer.getFeatures(QgsFeatureRequest().setFilterExpression("ID LIKE '{0}-%'".format(rang[i])))
    litp = list(itp)
    print 'nombre de pieds dans le rang ' + rang + ' : ' + str(len(litp))

I runned it and it works for first loop, first part of second loop and then it crash :

nombre de rang : 60
rang : 1
nombre de pieds dans le rang 1 : 6
rang : 2

exec((script), ns)
          File "<string>", line 34, in <module>
        IndexError: string index out of range

How do I fix this?


texts:
I have a colleague producing maps with the following Python Basemap projection:

from mpl_toolkits.basemap import Basemap
m = Basemap(width=6000000,height=3900000,
        resolution='l',projection='eqdc',\
        lat_1=45.,lat_2=55,lat_0=38,lon_0=-97.)
print(m.proj4string)

With a proj4 string of:

+y_0=1949999.9999999986 +units=m +lon_0=-97.0 +x_0=3000000.0000000014 +proj=eqdc +lat_0=38.0 +R=6370997.0 +lat_2=55.0 +lat_1=45.0

I'd like to match this in GRASS.

 g.proj -c proj4="proj4string from above" location=EQDC

Gives an error:

WARNING: Datum <unknown> not recognised by GRASS and no parameters found
Location <EQDC> created

I'm used to working in GRASS and using projections defined by GeoTIFF files, reprojecting, etc. I've created custom views in Python like above. But have not yet figured out how to tell GRASS to create something with the detail that I can in Python.


texts: How to use layers from trusted servers in ArcGIS Online?

texts:
We have a service on a secure server. It's available through a proxy; the service layers from the proxy can be added to an ArcGIS Online map and viewed, but not edited. It throws an error message saying they are not edittable because they are not available for arcgis.com. It's the same when I try to use the service directly - I get prompted for password, but I still can't edit.

Documentation seems to suggest trusting the service to bypass the proxy. However, all we have achieved yet was to avoid being prompted for password - however, there's no authentication working, so we can't get to the service neither. What are we missing?

The service supports CORS and allows our arcgis.com account (Access-Control-Allow-Origin: organization.maps.arcgis.com). After logging into the original server through REST API, I could add few features, so the REST endpoint is fully working.

EDIT: we use basic (HTTP Basic Auth or Web-Tier authentication), not token-based authentication. Sorry for not mentioning it in the beginning.


texts: Error in Google Earth Engine:Error generating chart: toFixed() digits argument must be between 0 and 100

texts:
I'm trying to do calculate land use's accuracy and kappa, this error happen in my code. I have tried a few things to sort this out，but I can't.
This is my code address：https://code.earthengine.google.com/?scriptPath=users%2Fjzhjk1067123%2Fyulin%3A%E6%A6%86%E6%9E%97%E5%B8%822015


texts: Changing central meridian in QGIS causes problems?

texts:
I want to change the central meridian of my world map (WGS84) in QGIS from 0 to -80. The problem is that this distorts the map every time I try it.

This is what I do:

go to Settings/Custom CRS and then click on Add new CRS
copy existing CRS (choosing WGS84 / World Mercator)
Edit its parameters (changing only the central meridian) to: +proj=merc +lon_0=-80 +k=1 +x_0=0 +y_0=0 +datum=WGS84   +units=m +no_defs
go to Project Properties/CRS and choose my custom CRS

And then the whole map is distorted. It works fine when I do it in ArcGIS.


texts:
As I have too much data in one shapefile (it is too long to load), I created from this shapefile 4 layers showing the level of environmental risk. However I would like these shapefiles to be subsetted to the original layer. So when I click on the + close to the original layer, the subsetted shapefiles appear.

How can I do that?


texts: Possible to drop NA geomType in readOGR?

texts:
I have a gdb file which is giving me trouble.
f = '/path/to/folder.gdb'
layers = ogrListLayers(f)
geom_list = lapply(layers, readOGR, dsn = f)

Gives an error on one of the layers:

Source: "f", layer: "troubleLayer"
Error in ogrInfo(...): Multiple incompatible geometries: wkbPolygon: 110284, NA: 493

It appears that some of the contained polygons are corrupted... I thought the require_geomType argument would fix this but alas:
readOGR(f, 'troubleLayer', require_geomType = 'wkbPolygon')

Gives a completely unexpected output:

Warning message in readOGR(...) : "eType not chosen"
[ repeat this message hundreds of times]
NOTE: keeping only NA wkbPolygon of 110777 features
Error: length(srl) > 0 is not TRUE

Is it possible to simply drop the apparent NA geometries and load this gdb into R?


texts: python script to convert shapefiles in a folder to .kmz's with the SAME shapefile name

texts:
Pretty new at working with python and am Stuck...I have a "for" Loop that looks in a folder and, based on a criteria (Homes valued over 500k), creates .kmz's with symbology for each of the .shp files in the folder. The script effectively creates the layers based on criteria and  applies symbology.  However I need help carrying the Name (fNAME) into the kmz's.  Can anyone help me complete this script?

Regards

Script:
~~~~~~~~~~~~~~~~~

import arcpy

from arcpy import env

# Variables
env.workspace = r"C:\Work\TEST\HOMES"

fcs = arcpy.ListFeatureClasses("D*", "Point")

HOMES_Lyr = r"C:\Work\TEST\HOMES\LAYER\HOMES.lyr"

arcpy.env.overwriteOutput=True

for fc in fcs:

    fNAME = arcpy.Describe(fc).basename+"_Green"
    arcpy.MakeFeatureLayer_management(fc,fNAME,"\"PRICE\" > 500000","","")
    arcpy.arcpy.ApplySymbologyFromLayer_management(fNAME,HOMES_Lyr)
    arcpy.LayerToKML_conversion(fNAME, ????.kmz, "12000", "false", "DEFAULT", "1024", "96", "CLAMPED_TO_GROUND")

~~~~~~~~~~~~~


texts: Calling geoprocessing service from another GP service arcpy.Importtoolbox?

texts:
I have a python script that calls an ArcGIS Server geoprocessing service toolbox with arcpy.Importtoolbox.  The GP service is on several separate machines behind a netscaler VIP.  Each machine is it's own environment, the machines are not a part of a site.  Now, when I run the python script from my local machine it runs just fine with no problems.  I can publish the script as a geoprocessing service, when I try to run this process from the rest endpoint, I get an error telling me it cannot find the tool:

  IOError: The toolbox file https://{netscalerURL}/arcgis/services;{folder}/{GP tool} was not found.

What I'm confused about is why the same script would work on my local machine but not as a geoprocessing service.  This is not the first time I've published a script like this but it is the first time we have tried to utilize a netscaler url as opposed to a specific machine in my call.

I think I have been able to identify the problem but I still don't have a fix.  The AGS machines are Linux based and run the GP service inside of Wine.  Wine does not recognize the SSL certifications on the machine which results in an invalid cert and blocked access.

Is there a way to incorporate the certificate into the arcpy.importtoolbox process using the ssl module?