texts:
Recently, I was looking for way to skeletonize buffer and besides solution with PostGIS and Postegre solutions there is algortihm v.voronoi.skeleton. I have tried to use it with different road/rail buffers as described in this tutorial (unfortunatelly in polish, though the best source I got).

Despite following instructions with different parametres and input layers I always get error with output and no layer is created.

2019-07-01T20:50:58     WARNING    Duplicate output OUTPUT_HTML_FILE registered for alg build_generic_plot
2019-07-01T20:50:58     WARNING    Duplicate output OUTPUT_JSON_FILE registered for alg build_generic_plot


texts:
Currently, I am having some problems with ArcGIS Online custom Geoprocessing tools.

I made a custom tool using ModelBuilder and published it on my portal. Even though they worked just fine with no problem whatsoever in ArcGIS Desktop, they seem to give out an error when running in Online.

The error is nothing more than "failed" and i was wondering if anybody knows how to get more details of what went wrong so I can try to fix I?

Currently I have nothing to go on.


texts: r.series GRASS algorithm in QGIS 3.2 output no-data image

texts:
I'm trying without success to calculate a standard deviation image from 4 input images. I already used the r.serie algorithm in 2.18 to calculate an averaged image. Know after upgrading to 3.2.2 on ubuntu 18.04 r.series seems not to work properly. the output is an image with value from 1.79769e+308 to -1.79769e+308.

Any ideas? The log is:

Cleaning up temporary files...

access: No such file or directory

ERROR: LOCATION </home/gobel/newLocation> not available

Starting GRASS GIS...

Executing </tmp/processing_c22010c5e1ed49c0ae5fc8fee008fb9d/grassdata/grass_batch_job.sh> ...

Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default

Projection information updated

Over-riding projection check

Reading band 1 of 1...

r.external complete. Link to raster map <rast_5b7ee239a136610> created.

Over-riding projection check

Reading band 1 of 1...

r.external complete. Link to raster map <rast_5b7ee239a139e11> created.

0..3..6..9..12..15..18..21..24..27..30..33..36..39..42..45..48..51..54..57..60..63..66..69..72..75..78..81..84..87..90..93..96..99..100

Checking GDAL data type and nodata value...

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

Using GDAL data type <Float64>

Input raster map contains cells with NULL-value (no-data). The value -nan will be used to represent no-data values in the input map. You can specify a nodata value with the nodata option.

Exporting raster data to GTiff format...

ERROR 6: SetColorTable() only supported for Byte or UInt16 bands in TIFF format.

2..5..8..11..14..17..20..23..26..29..32..35..38..41..44..47..50..53..56..59..62..65..68..71..74..77..80..83..86..89..92..95..98..100

r.out.gdal complete. File </tmp/processing_c22010c5e1ed49c0ae5fc8fee008fb9d/9331a0aec53041829431d0b6eb7c1af7/output.tif> created.

Execution of </tmp/processing_c22010c5e1ed49c0ae5fc8fee008fb9d/grassdata/grass_batch_job.sh> finished.

Cleaning up temporary files...

Execution completed in 24.22 seconds
Results:
{'output': <QgsProcessingOutputLayerDefinition {'sink':/tmp/processing_c22010c5e1ed49c0ae5fc8fee008fb9d/9331a0aec53041829431d0b6eb7c1af7/output.tif, 'createOptions': {'fileEncoding': 'System'}}>}

Loading resulting layers
Algorithm 'r.series' finished


texts: connecting Arc 10.3 to a PostGIS database "the coordinate system identifier is invalid"

texts:
Has anyone had any success directly connecting Arc 10.3 to a PostGIS/postgres database, where the geometry type is 'geometry'?

I am able to connect and view the tables, but I am unable to view them - I receive the error from ArcCatalog "the coordinate system identifier is invalid".

I have created a geodatabase out of the data using 10.3 and a provisioning file to create sde.db_tune table, but that has made no difference (FYI, geometry_storage in db_tune is "PG_GEOMETRY"). Once I created the SDE tables, I attempted to connect using 10.2 but also received the same error message.

SRID is set per QGIS table properties (QGIS views the data perfectly, by the way).

I'm using PostGres 9.3 and PostGIS 2.1. I'm using extensions "postgis" and "postgis_topology" in my public schema, so minimal table massaging appears to be necessary (like manually updating geometry_columns).

Any ideas?


texts: Displaying polylines in different color / Color gradient

texts:
Though i am using QGIS, yet I am not confined by tool (can use an entirely different tool or a plugin for QGIS) usage.. Only aim is to get the following:

I am reading various shapefiles in QGIS and I wish to have different colors in the same shapefile (or layer) for -

Each polyline
Every n number of polylines
A chunk of n points in a single large polyline.


texts: Trying to project GeoPandas geodataframe geometry to new CRS is not working as expected

texts:
I am working with GeoPandas and a Point GeoDataFrame. I have data points in WGS84 (EPSG=4326) that are located in New Zealand. I am trying to project them to a local UTM, but I cannot get the geometry of the GeoDataFrame to update with the new UTM coordinates.
Here is an example GeoDataFrame:
import geopandas as gpd
from pyproj import CRS
long = [172.083916, 172.086549, 172.093403]
lat = [-43.594391, -43.594758, -43.595249]
gdf_pt = gpd.GeoDataFrame(geometry=gpd.points_from_xy(long, lat))
gdf_pt.crs = 4326
display(gdf_pt)
print(gdf_pt.crs)

Failed Attempt #1 to project to EPSG=2193
myproject = gdf_pt.copy()
myproject.crs = 2193
print(myproject.crs)
display(myproject)

Failed Attempt #2
myproject['geometry'] = myproject['geometry'].to_crs(epsg=2193)
display(myproject)

I also tried it with EPSG=32759, another suitable NZ UTM... same issue
What am I doing wrong?


texts:
I have a field team using a Trimble Geo XH from 2008.

My usual workflow is to have them email me the raw GPS files, which I copy into a "GIS Folder" for import to Pathfinder where I process it.  I then transfer the .cor file to that GIS folder.  The output is 8 files that together make the GPS data file. I email the those 8 files to the team to have them copy/paste them on the the GPS unit.

They pointed out to me that the files I sent do not contain a ".CAR" file (no, not ".cor") but there is one currently on the GPS unit. I do not believe I have the ability to recreate a CAR file

What information is in the ".CAR" file?  Can we remove it? Is it important?  Is it optional like .prj .sbn  and other components of a shapefile are?

Edit: I had the team remove the .CAR file from the GPS before copy/pasting the 8 GPS files over. Data collection appeared to go fine all day. By the end of the day when they  backed up the data, the .CAR file had returned, however H-Star processing options are no longer available!


texts: Failing to use the correct Con expression in raster calculator in ArcMap

texts:
So, I am trying to use the Con statement in ArcMap to add a certain value to a raster dataset (raster2) if the data is the same as an other dataset (raster1), if not I should just stay the same. I tried to use the expression:

Con(("raster1" == 1), ("raster2" + 1.5), "raster2")

But it seemed to null out any values that where not within raster 1, and therefore loses a lot of values and data.

Does anyone have a solution?


texts: How to retrieve data from a related table by using ArcGIS JavaScript API 4.1?

texts:
I have a feature layer and a related table. I want retrieve information from the related table and show them.

Which API is referenced when using ArcGIS JavaScript API 4.1.?


texts: Displaying "COUNT" from a point-shapefile in ArcGIS 10.2.2

texts:
I've tried a number of ways, labeling, unique values, etc.  I know I'm right there but can't get it, it will just display features by the "NAME" field.  Also if I try to use symbolize by quantities, only fields I get are numerical I'm guessing, but cannot be aggregated to a total there.  I cannot find the Field Calculator either, which I thought I could maybe add a field and have a "SUM" field added and then label it that way.

End Goal:

I have point features in various states, and the application needs to display the number of those features within each state.


texts:
I am using Leaflet.Draw to draw a circle.

var circle = L.circle(new L.LatLng(latitude, longitude), number, circle_options).addTo(featureGroup);

I need the coordinates of the circle as a polygon (LatLng[]).

Is it possible?


texts: Assigning OpenLayers Layers to different layer trees 'automatically'

texts:
I am developing an OpenLayers 3 web application were I would like to incorporate multiple layer trees. Instead of manually adding each layer to a layer tree I would like to do it using the method below.

I have this JSFiddle set up that it has a map with multiple layers and one Layer Tree. I would like to duplicate the layer tree multiple (4) times and have a parameter (as per the one in each layer of the example) that then assigns a layer to different layer trees.

I imagine it involves an extra if statement when building the layer tree at the bottom of the code, similar to the slider option already in there, but I don't seem to be able to make it work.

The solution should also ideally remove the "all layers" group at the top of each Layer Tree.


texts: Convert all the polygons of a FeatureCollection to editable polygons in GEE

texts:
I need to convert all the polygons of a FeatureCollection in editable polygons (like the polygons you can draw with the drawing tools) so I can modify each of their shape if needed. The polygons represent burned areas during 1 year of a fire season.

var table = ee.FeatureCollection("users/cecilianaval/Poligonos_SC_BAMS_2018");
print(table);

//Center and display the table, zoom = 8
Map.centerObject(table, 8)
Map.addLayer(table, {}, 'BAMS_2018');

Code Editor script


texts: R use data in SpatialGrid to create SpatialPoints

texts:
I am trying to make a grid of points to overlay on a raster of vegetation of Kansas, U.S.A. My idea is to use the points as a sort of centroid for the raster pixels, which are 30m resolution. When I have the points, or "centroids", I want to use getValues() to capture some attributes of the raster and add the values from getValues() to the points. The steps I have taken so far, from Bivand et al.:
> bbox(ks08)
         min       max
s1  229687.8  889807.8
s2 4094357.4 4434227.4
> cs <- c(30, 30)
> cc <- bbox(ks08)[,1]+ (cs/2)
> cd <- ceiling(diff(t(bbox(ks08)))/cs)
> p4k <- CRS(proj4string(ks08))
> p4k
CRS arguments:
 +proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
+no_defs
> ks_grid <- GridTopology(cellcentre.offset = cc, cellsize = cs, cells.dim = cd)
> ks_sg <- SpatialGrid(ks_grid, proj4string = p4k)
> summary(ks_sg)
Object of class SpatialGrid
Coordinates:
         min       max
s1  229687.8  889807.8
s2 4094357.4 4434227.4
Is projected: TRUE
proj4string :
[+proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
Grid attributes:
   cellcentre.offset cellsize cells.dim
s1          229702.8       30     22004
s2         4094372.4       30     11329
> str(ks_sg)
Formal class 'SpatialGrid' [package "sp"] with 3 slots
  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
  .. .. ..@ cellcentre.offset: Named num [1:2] 229703 4094372
  .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
  .. .. ..@ cellsize         : num [1:2] 30 30
  .. .. ..@ cells.dim        : int [1:2] 22004 11329
  ..@ bbox       : num [1:2, 1:2] 229688 4094357 889808 4434227
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "s1" "s2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=utm +zone=14 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

When I view ks_sg, part of the output include SpatialPoints and then the output chokes the memory (249,000,000+ points) for a few moments while it prints all the points' coordinates. I don't know what to do next to get those points into a vector format, like SpatialPoints, that would be the equivalent of a ArcGIS shapefile. Making the vector of points is my goal at the moment.
After the reply from Dave-Evans: I did try rasterToPoints before this, but ran into memory issues. Setting up the grid seemed like my best alternative.
After reply from Jeffrey Evans: My ultimate goal is to set up a very large data frame/data.table that has a column for each year of cover. The column will be vegetation classes for the year for every pixel, like a cbind(getValues(ks08),getValues(ks09),...). Each row is a pixel and each column is a year and the cell for a pixel and year would have the vegetation cover at that location and time.


texts: Settings of Address Locator within ModelBuilder or ArcPy?

texts:
I am creating a model that processes a large road dataset, then creates an address locator from the roads. The Create Address Locator tool does not allow me to set options like Spelling Sensitivity or Minimum Candidate Score.

Is it possible to set these options in ModelBuilder (or at least, within an ArcPy Python script)?

I'm using ArcGIS 10, but I'm also interested in whether it's possible in 10.1.


texts: Extracting Spatial Points that match Raster Value using R?

texts:
I have a question concerning the extraction of spatial points that match a certain value of an additional raster in R. The operation itself should actually be pretty simple, but I'm stuck somehow. (I am using R at a "lower intermediate" level).

I am using two data sets. The first one is a shapefile containing point data.

training1 <- shapefile ("re_27042015projcrop_all.shp")

The second one is a "binary mask" raster object consisting solely of 0 and 1 values.

radar1 <- brick ("binary_01_180115_240415.envi")

How do I extract all points whose coordinates match raster pixels with a value of 1?

I tried some logical expressions like

output <- training1@coords & radar1==1

and various raster functions like crop, mask etc., but none seemed to work.


texts: Plotting raster and vector data on same Plot

texts:
I have a raster and croped it by vector boundary. Vector reprojected to crs of the raster, but when I plot it does not show on the same figure?

Something with plotting extend, but how to make them match even if the coordinate system the same?

# Process

# Defines path to data
fire_boundary_path = os.path.join("...200_dd83.shp")
naip_2017_path = os.path.join("...70902.tif")

fire_boundary = gpd.read_file(fire_boundary_path)

# Reproject fire boundary in the CRS as NAIP
fire_boundary_pr = fire_boundary.to_crs(epsg=26913)
fire_boundary_pr.crs

# Open/crop the rasters data sets
with rio.open(naip_2017_path) as naip_2017_src:
    naip_2017_crop, naip_2017_crop_meta = es.crop_image(naip_2017_src, fire_boundary_pr)

# To Plot
fig, (ax) = plt.subplots(figsize=(14, 14))
ep.plot_rgb(naip_2017_crop,
            ax=ax,
            rgb=[3, 0, 1],
            extent=naip_2017_crop_meta,
            title="Hw Plot#: NAIP 2017 Post Fire CIR Image")

# Including the Cold Springs fire boundary
fire_boundary_pr.plot(ax=ax, color="none",  edgecolor="blue", linewidth=3.1)

#plt.show()


texts: Want to use a map widget with Ozone Widget Framework (OWF)

texts:
I have written some widgets for OWF.  They work well and I can easily send and receive intents and I can send messages on channels.  That's great.

I can also send Intents and messages for the CMAPI (Common Map API) which is supposed to be the standard way of saying, "please draw this on a map".

The only thing I can't do is find an implementation of CMAPI.  It's wonderful that there is a standard for interacting with maps, but has anyone tried to implement it before?

I did find the ESRI OWF widgets called Esri-Ozone-Map-Widget on github.  That's cool and sounds like exactly what I want, except it says:

  CMWAPI map.feature.plot messages are not supported by the map widget

What good is a map widget if you can't plot features on it?  I have no idea.  If I were implementing a map widget, the first thing I would work on is map.feature.plot.  But I'm not a GIS expert.  Maybe there are many things people want to do with maps other than plotting features on them.

I also found Boeing Meridian project.  That works great and allows me to open a map as a widget.  It claims that it responds to CMAPI intents but I was never able to get that working.

Is there anything else?  Surely someone has thought, "hey, let's try using maps to display geospatial data within OWF"?  I can't be the first one to consider this idea...

There is a simple Google Maps widget that comes with it as a demo, but it only supports plot.  I can't even remove a marker after I have plotted it.  It's not for any serious use.


texts: Want to insert a new field to CSV with filename as attribute

texts:
I have thousands of CSV files with the approximately the following structure:

csv_20160601T004301Z:

name,x,y,
name1,23,66
name2,23,67

csv_20160601T005001Z:

name,x,y,
name1,23,66
name2,23,67

etc.

As you can see, the filename contains a timestamp. I'd like to add a new field to each of the csv files and populate the field with the filename, like so:

csv_20160601T004301Z:

name,x,y,timestamp
name1,23,66,csv_20160601T004301Z
name2,23,67,csv_20160601T004301Z

csv_20160601T005001Z:

name,x,y,timestamp,
name1,23,66,csv_20160601T005001Z
name2,23,67,csv_20160601T005001Z

etc.

Of course, the "csv_" prefix is not really needed, just wanted to illustrate the idea. Any thoughts on how to do this? Is this something that can be achieved with Notepad++ or similar? Or will this require a more sophisticated approach?

I have very, very little experience in coding/programming, so I may need the potential ideas explained in the simplest way possible.

edit: SOLVED. I found a solution that helped me. I used the following in cmd:
for /f %a in ('dir /b .csv') do for /f "tokens=" %b in (%a) do echo %b,%a >> all.csv

Source: https://community.spiceworks.com/topic/322603-merge-csv-files-and-create-a-column-with-the-filename-from-the-original-file


texts:
I've created a leaflet map using maptiler.

It's created a leaflet HTML file, which links to PNG images on my desktop and displays the image.

How do I embed the map onto a blog such as tumblr?


texts: Does anyone know where I can find Stream (rivers, waterbodies) codes for the State of Utah?


texts: error while running CGAl pacakge examples

texts:
I'm try to run "map_3_with_colored_facets" example which included in CGAL package
but I faced some errors , it cannot read "Combinatorial_map.h" methods

-cm.make_combinatorial_hexahedron() isn't a member
-  cm.attribute<2> isn't a member
- cm.info<2>  isn't a member ....

#include <CGAL/Combinatorial_map.h>
#include <CGAL/Combinatorial_map_constructors.h>
#include <CGAL/Cell_attribute.h>
#include <iostream>
#include <algorithm>
#include <cstdlib>
struct Sum_functor
{
template<class Cell_attribute>
void operator()(Cell_attribute& ca1,Cell_attribute& ca2)
{ ca1.info()=ca1.info()+ca2.info(); }
};

struct Divide_by_two_functor
{
template<class Cell_attribute>
void operator()(Cell_attribute& ca1,Cell_attribute& ca2)
{
ca1.info()=(ca1.info()/2);
ca2.info()=(ca1.info());
}
};
struct Myitem
{
template<class CMap>
struct Dart_wrapper
{
typedef CGAL::Dart<3, CMap> Dart;
typedef CGAL::Cell_attribute<CMap, int, CGAL::Tag_true,
 Sum_functor, Divide_by_two_functor> Facet_attribute;
 typedef CGAL::cpp11::tuple<void,void,Facet_attribute> Attributes;
 };
};
typedef CGAL::Combinatorial_map<3,Myitem> CMap_3;
typedef CMap_3::Dart_handle               Dart_handle;
int main()
{
CMap_3 cm;
// Create 2 hexahedra.
Dart_handle dh1 = cm.make_combinatorial_hexahedron();
Dart_handle dh2 = cm.make_combinatorial_hexahedron();
// 1) Create all 2-attributes and associated them to darts.
for (CMap_3::Dart_range::iterator
   it=cm.darts().begin(), itend=cm.darts().end();
   it!=itend; ++it)
 {
 if ( cm.attribute<2>(it)==NULL )
  cm.set_attribute<2>(it, cm.create_attribute<2>());
 }
 // 2) Set the color of all facets of the first hexahedron to 7.
 for (CMap_3::One_dart_per_incident_cell_range<2, 3>::iterator
   it=cm.one_dart_per_incident_cell<2,3>(dh1).begin(),
   itend=cm.one_dart_per_incident_cell<2,3>(dh1).end(); it!=itend; ++it)
 { cm.info<2>(it)=7; }
 // 3) Set the color of all facets of the second hexahedron to 13.
 for (CMap_3::One_dart_per_incident_cell_range<2, 3>::iterator it=
   cm.one_dart_per_incident_cell<2,3>(dh2).begin(),
   itend=cm.one_dart_per_incident_cell<2,3>(dh2).end(); it!=itend; ++it)
 { cm.info<2>(it)=13; }
 // 4) 3-Sew the two hexahedra along one facet.
 cm.sew<3>(dh1, dh2);
 // 5) Display all the values of 2-attributes.
 for (CMap_3::Attribute_range<2>::type::iterator
   it=cm.attributes<2>().begin(), itend=cm.attributes<2>().end();
   it!=itend; ++it)
  {
   std::cout<<cm.info_of_attribute<2>(it)<<"; ";
 }
std::cout<<std::endl;
// 6) Insert a vertex in the facet between the two hexahedra.
cm.insert_cell_0_in_cell_2(dh2);
// 7) Display all the values of 2-attributes.
for (CMap_3::Attribute_range<2>::type::iterator
   it=cm.attributes<2>().begin(), itend=cm.attributes<2>().end();
   it!=itend; ++it)
{
  std::cout<<cm.info_of_attribute<2>(it)<<"; ";
}
std::cout<<std::endl;
cm.display_characteristics(std::cout);
std::cout<<", valid="<<cm.is_valid()<<std::endl;
return EXIT_SUCCESS;
}


texts: SQL expression for largest number in Table

texts:
With "select by attributes" I would like to select the largest number in a column of an attribute table. How would be the SQL expression for that?

I don't want to select it by hand, because I will use it later in a python script (arcpy.SelectLayerByAttribute_management()).

I am using ArcGIS 10.4 and the files are within a Personal-Geodatabase.


texts: Want to access android location services in a leaflet.js project running in a webview

texts: I am running leaflet.js project on my android app using a webview. Everything is great but for current user location leaflet uses GPS because it runs in a webview and as we know webview is extention of browser on android. Now GPS is not very reliable. I want to access android sensors of accelerometer or geomagnetic sensor to determine user position not GPS. Is there a way to achieve this?


texts: Using the ArcGIS 'select by location' tool in ModelBuilder?

texts:
I have a model that outputs a shapefile which shows the extent of flooding. I also have a layer of buildings. I want to see which buildings are affected by different flood scenarios.

When I use Select by Location, my results look like this, which is what I want:

However, Select by Location isn't something that can be added to a model, so I'm looking for a workaround. Unfortunately, using the Clip and/or Intersection tools, I get many small artifact polygons. I was wondering how to add Select by Location to a model, or if there was some way to perform an overlap/within type of intersect rather than the default.


texts:
I just can't get the latest GDAL release running with Visual Studio 2019.
Creating the VS project is no problem. Pretty straight forward. But when I try to build the project in VS in fails and gives me this error:
PROJ_INCLUDE should be defined. PROJ >= 6 is a required dependency

Of course I found this question here: https://stackoverflow.com/questions/57600419/error-building-gdal-proj-include-should-be-defined-proj-6-is-a-required-dep
But that doesn't solve the problem. I tried downloading proj with OSGeo4w and I modified the nmake.opt file accordingly.
I tried every link I could find on google but with no luck. The GDAL community seems to be aware of that but maybe just don't care?
What can I do? Or do I have to use an older GDAL 2.x version (which seems to work) ...


texts: How can I use regex in the AttributeCreator Advanced Editor in FME 2014?

texts:
I want to use regex to simply distinquish between POLYGON and MULTIPOLYGON using the regex term "\bPOLYGON"

Apparently I can use the Advanced Editor and use

@FindRegEx(<STRING>,<STRING>)POLYGON

But what do I enter here for STRING? In the documenation it states

FindRegEx(string str, string regExp, [int startIdx], [bool caseSensitive], [-1])

So I guess "string regExp" is the expression "\bPOLYGON", and I would guess "string str" means the source text....but this is already defined as the "left Value".

Can anyone give me an example of how to use this?


texts: Placing random line segments within polygon?

texts:
I would like to generate multiple line segments of a fixed length, a random bearing, and both segment endpoints within the polygon.
It's not clear from the example, but it's fine if the segments intersect each other or span the indented features as long as both endpoints are in the polygon.

I also want the attributes table to include lat/lon for each segment endpoint and the azimuth between the endpoints.

Is this possible?


texts: Spatial Adjustment of shapefile not possible

texts:
I am trying to spatially adjust a shapefile containing several polygons to another shapefiles where spatial data is correct. Both should share exactly the same extent. However, when I am using spatial adjustment it never lets me accurately match the two shapefiles based on the displacement links. No matter how many links I set. An obvious shift always remains.

The two shapefiles can be seen HERE:
https://drive.google.com/drive/folders/10mcprmxrLFzUnYRYILZwE5hAO6bcuNGe?usp=sharing

"PopDensity" is the one that needs to match the correct "statistische_gebiete"-shapefile.

Any clue what I am doing wrong? Or is my understanding wrong that the spatial adjustment scales, rotates and shifts whole shapefiles?


texts: Seeking symbology (.lyr?) for USGS geologic maps?

texts:
I work a lot with USGS geologic maps and always have the same problem. I am able to download the files from the website (in .e00 format), and make a nice shapefiles out of them. My problem is that there are no .lyr files with symbology, and assigning colors to polygons manually takes ages. The USGS have always great databases and metadata, so it is hard to believe that they wouldn't attach any file which would help with coloring polygons... Does anybody know by chance if there is any symbology "hidden" in files like .e00, .ps, .aml, .met? The database goes with those and I am not sure what are they for...

My data:
http://geo-nsdi.er.usgs.gov/metadata/map-mf/2337/metadata.faq.html

The map looks like this:
http://pubs.usgs.gov/mf/2000/2337/mf2337f.pdf


texts: How to put a Google Map image (or any other map layers) under my UTM cordinates plot uding matplotlib and Python?


texts: Command line error and missing dependency in QGIS

texts:
I have two problems with QGIS 2.6.1

1) When i am trying to open command line i get error:

An error occurred during the execution of Python code:

Traceback (most recent call last):
File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\ProcessingPlugin.py", line 134, in openCommander
    self.iface.mapCanvas())
  File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\gui\CommanderWindow.py", line 48, in init
    self.commands = imp.load_source('commands', self.commandsFile())
UnicodeEncodeError: 'ascii' codec can't encode character u'\xf3' in position 17: ordinal not in range(128)

2) When i am trying tu run external algorithm SAGA, GRASS whatever i get message:

Missing dependency.This algorithm cannot be run :-(
This algorithm requires SAGA to be run.Unfortunately, it seems that SAGA is not installed in your system, or it is not correctly configured to be used from QGIS
Click here to know more about how to install and configure SAGA to be used with QGIS


texts:
I have a map with a "strange" projection (it's basically hand-painted, not an exact one). I can easily georeference it in QGIS using the Georeferencer plugin and the Thin Plate Spline transformation with a set of some 2000-3000 control points into simple WGS 84 (EPSG:4326) "latlong" projection and process/enhance it further with other data, like from Natural Earth or OpenStreetMap.

Now I'd like to perform the inverse operation: Transform one of those raster layers (or rasterised vector layers) back into the hand-painted quasi-projection, ideally using the same set of control points, so that the images match when overlaid in a graphic editor. How do I do it?


texts:
I use gdaldem to calculate the slope of a given region. This is the command I use:

gdaldem slope -s 111120 -p -compute_edges  N34E069.hgt N34E069_slope.tif

(I've got the SRTM data from http://viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm)

Some values in the raster band however have the noData value -9999. In the documentation I just came accross this line:

  The value -9999 is used as the output nodata value.

But, what is the reason for getting a nodata value? When I create a color-relief map with the following color code:

-9999 255 0 0
0 255 255 255
90 0 0 0

I get the following output:

It seems like the noData values occur on the peaks of mountains. But on the left in the flat area I also get some noData values.

How can I fix this problem?


texts: Merging vector with many polygons to single vector polygon in QGIS?

texts:
I am trying to add the boundaries of one shapefile (A) that are contained within a single polygon to the shapefile (B) of that single polygon.  Both are at the same projection.  I have tried clipping, intersecting, union, etc. but have not been able to accomplish it.

example http://radslav.slavniklabs.net/dwnld/problem.jpg

Can anyone provide step-by-step instructions on how to accomplish this using QGIS 2.10.1-Pisa?

When I use Geoprocessing/Clip A(DouglasCountyNE2) with B(sitearea2) the dialog looks like this:

dialog http://radslav.slavniklabs.net/dwnld/clipdialog.png

The result is testomaha.shp but the attribute table is empty:

results http://radslav.slavniklabs.net/dwnld/clipresult-nothing.png

If I use Geoprocessing/Difference the result includes areas beyond the area B as seen below:

difference http://radslav.slavniklabs.net/dwnld/difference.png

This is because it did not clip the areas in A to match B.

Here is a link to a zipfile with both shapefiles in it:

http://tinyurl.com/qdf7dhu


texts: Reconstituting features from E00 format, without ArcGIS

texts:
Recently I used OGR to convert an ESRI E00 format file to shapefiles. I am dealing with vector data, many polygon features per file. I was dismayed to find that none of the resulting shapefiles(you get one shapefile for each sub-file within the E00 file) 'reconstitute' the original features. According to the only E00 format documentation on planet earth there is an ARC file which describes all the arcs used for all features, and a PAL file which lists which arcs are used to describe each feature, and how the arcs should be combined.

As far as I can tell, the only way to get the actual features out of a E00 file is to use ArcGIS which can import the E00 and export all the shapefiles within. Is there any other open source possibility that will automatically extract the features? Or is the only solution to convert the E00 to shapefiles, and then write some code to brute force the features from the resulting ARC and PAL shapefiles i.e. parsing the ArcIds attribute from the PAL and building up a feature from the listed arcs?

End of main question.

Note that:

ogrinfo -al

and looking at the resulting ArcIds in the PAL section yields such gems as:

ArcIds (IntegerList) = (12:-58,164,197,220,275,274,-319,-344,-346,-329,-265,-203)
ArcIds (IntegerList) = (11:216,267,291,-292,257,0,-237,0,-239,0,-158)

where 0 indicates the following arc describes a hole, and a negative index indicates that the arc must be reversed. Not sure how 12:-58 would be interpreted. I assume

12,11,...,1,-1,-2,...-58

The format also indirectly implies that all holes must be comprised of a single closed arc.


texts: Use define fixed size rotated SVG symbol controlled by Atlas

texts:
I'm working with QGIS 3.10 under FreeBSD 12.1 and trying to use a North direction SVG arrow symbol which is controlled by an atlas. The atlas is driven by a set of map sheets arranged along a winding corridor.

Each sheet has a fixed size and scale as well as the map is turned into a axis parallel chart by the sheet (atlas driven). For each sheet I defined a North arrow rotation parameter corresponding to map rotation ("Composer -> Item Properties -> North Arrow -> Image rotation").

Unfortunately the size of the North arrow is automatically adopted for each rotation due to the corresponding rotated angular dependent bounding box ("Main Properties -> Resize Mode").

Example North Direction

Example North-East Direction

If I use of main "Properties -> Clip" (which I guess is the right parameter for a circular symbol which is centered an cannot be clipped), the rotation by the atlas is turned off.

Which parameter do I've to use get a rotated equal sized atlas controlled symbol under the "Composer -> Item Properties -> North Arrow" dialog.


texts: Possible to use same database for mapnik and pgRouting?

texts:
I am already running an OGCServer (mapnik) with an postgis database.

I want to also run pgRouting on preferably the same database, so that I do not have to maintain two separate databases.

I understand that pgRouting adds/changes some requirement to the data model. My question is basically - does it change the requirements in an incompatible way, or does it just add requirement.

(As you might guess; i'm still quite a newbie at pgRouting)


texts: Within a GeoServer framework, does anybody know how I might best create/cache tiles based on layers stored in PostGIS?  I am aware of Mapnik and have created a simple TMS using python (+ Django) but this seems 'clunky' and I'd like to know if there are any better ways to go about it.


texts: How do I get GDAL_translate to output better resolution

texts:
I have a set of GeoPDFs that I converted to GeoTiffs. However, the end results are horrible. They are very pixelated. The command I passed was

gdal_translate -of GTiff inFile outFile -co COMPRESS=NONE -co JPEG_QUALITY= 100.

Is there a setting that I am missing because some of the outputs are unreadable. By comparison the outputs I get from Global Mapper 15 are super sharp. Sadly there is a whole other world of problems with GLmapper though.

I was interested in the resample settings but figured that the default method of nearest neighbor would be the sharpest.


texts: Calculating meshes average within boundary?

texts:
I have mesh data of elevations. Each mesh has a value of elevation (e.g. 87, 60, 33, ...).

And I also have data of administrative boundaries.

I want to calculate the average of each administrative area's elevation, but I don't know how to connect two shapefiles (mesh data and boundary data) because these two are different shape.

Could someone tell me how to?


texts: Rotating rasters using GDAL Geotransform (and ArcGIS Desktop)?

texts:
I am working on a GDAL driver for a proprietary grid-format. In some cases, the grids are rotated in order to get smaller file sizes. In my GDAL geotransform array, I do the following (showing a hardcoded 30 degrees angle):

rot_deg = 30;
rotation = rot_deg * deg2rad;
poDS->dGeoTransform[0] = originX;
poDS->dGeoTransform[1] = cos(rotation) * cellSizeX;
poDS->dGeoTransform[2] = -sin(rotation) * cellSizeX;
poDS->dGeoTransform[3] = originY;
poDS->dGeoTransform[4] = sin(rotation) * cellSizeY;
poDS->dGeoTransform[5] = cos(rotation) * cellSizeY;

where originX and originY is the center of the corner cell.

In addition, I got an identical file in TIFF format, with same min/max + cell values.I then apply an ArcGis geotransform to the .tiff, rotating it 30 degrees aswell.

When adding both files as layers, the tiff is 100% correctly aligned to the pivot point, while my GDAL format is slightly displaced.

The screenshot shows the overlapping grids, with only the corner pixel visible. The cyan point is the pivot point. The dark pixel on top is the (correct) tiff, while the grey pixel beneath is from my own driver. The misalignment is the same over the whole grid.

Anybody got an idea what could be wrong ?

In the ArcGis geotransform, it is possible to enter a pivot point, which is the case here (center of the corner cell). However, it seems like GDAL's geotransform is not expressive enough (?) to handle this, hence it uses the lower left corner.


texts: QGIS 3.0 is showing a point's coordinates in the original spatial reference system instead of the project's spatial reference system

texts: I obtained GPS coordinates of a point to within a couple centimeter accuracy using RTK GPS in the NAD83 Spatial Reference System. This information is included in the QGIS layer file metadata, and the properties properly shows the spatial reference system as NAD83. However, I have my project spatial reference system set to WGS84 spatial reference system, and when I hover over the point, the coordinates displayed are the NAD83 coordinates instead of the WGS84 coordinates (which are different by about 1 meter in north america.I am using the UTM16N projection for both spatial reference systems, so the coordinates are displayed in meters. How can I get the map to display the transformed coordinates in the project's spatial reference system so that I can compare the coordinates of points with different spatial reference systems? ArcGIS displays the coordinates using the spatial reference system of the "view", and it seems misleading for QGIS not to be doing this.


texts: Link between node->building->street in OSM

texts:
Data: Raw OSM data (.osm XML file).

Target: Connect between house number represented by node, to street name represented by way.

What I did:

In www.openstreetmap.org I can explore elements. So I point the node representing house number:

Then I searched by it's osm_id and found:

<node id="MY_NODE_ID" visible="true" version="1" changeset="..." timestamp="..." user="..." uid="..." lat="..." lon="...">
    <tag k="addr:housenumber" v="MY_HOUSE_NUMBER"/>
</node>

Then I searched this node and as expected I found it part of the building represented as way:

<way id="MY_BUILDING_ID" visible="true" version="1" changeset="..." timestamp="..." user="..." uid="...">
    <nd ref="MY_NODE_ID"/>
    <nd ref="..."/>
    <tag k="building" v="yes"/>
</way>

Finally I expect MY_BUILDING_ID to be part of another way representing the street. But MY_BUILDING_ID appears only once in the XML file. So I found street's osm_id and found it in the XML as:

 <way id="MY_STREET_ID" visible="true" version="10" changeset="..." timestamp="..." user="..." uid="...">
  <nd ref="..."/>
  <nd ref="..."/>
  <tag k="highway" v="residential"/>
  <tag k="name:en" v="MY_STREET_NAME"/>
  <tag k="source" v="user via ..."/>
 </way>

But none of the <nd ref="..."/> isn't MY_BUILDING_ID, and more over none of the <nd ref="..."/> appears again in the XML file.

Conclusion: How to link between known MY_NODE_ID to MY_STREET_ID so later it's easy to make the link between MY_HOUSE_NUMBER and MY_STREET_NAME?


texts: More precise control over QGIS symbology classification

texts:
In QGIS, I would like to classify my graduated symbology by standard déviations, as such:

. However, I find that the behaviour seems erratic or hard to control, and I can't choose the number of classes like I do with other classification methods (i.e. quantiles). Is there a way to have a more precise control over how many classes (and deviations from the mean) are used to classify the data set?

For instance, I would like to classify this distributio, by steps of 0.5 std devotions. So that the scale would range from < -1.5 Std Dev to > +3.5 Std Dev :


texts: Raster Analysis in Geoprocessing Service published from ArcMap

texts:
I have prepared a tool for Masking the specific portion of the raster dataset using the feature dataset. It works fine in ArcMap but shows the following error while running it on the server.

I have used shapefile for the feature dataset and Tiff file for the raster dataset.


texts:
QGIS does not have the coordinate reference system I need, So I decide to define a new custom Geographic coordinate system.To define a CRS, i have srs.db & .wkt file.
I place these file as
(a) Copy SRS.db in the folder and paste this file in C:\Program Files\QGIS 3.4\apps\qgis\resources
(b) Copy ESRI_EXTRA.wkt: paste this file in C:\Program Files\QGIS 18\share\gdal

After successful data import my custom CRS shown in predefined coordinate reference system:

Problem is when i try to these CRS for my project the CRS is Selectable. It shows check in No projection.


texts: How to calculate the slope of my polygons using DEM

texts: I downloaded elevation data (DEM) of my study area from NASA and I converted it to slope. How do I calculate the slope of each farm at my study area? I have a shapefile containing the polygons of the farms


texts: What does "foo" stand for in this PostGIS command?

texts:
The help page for the PostGIS function ST_PixelAsCentroids has this line of code as a usage example:
SELECT x, y, val, ST_AsText(geom) FROM (SELECT (ST_PixelAsCentroids(rast, 1)).* FROM dummy_rast WHERE rid = 2) foo;

What is the word "foo" standing in for in this example?


texts: load map parameters after opening new window

texts:
I am trying to open a new window from an existing window, which will contain a map with an extent derived from the original window.  In other words, the user will zoom to an area in the original window, click a button, then a new window will open with a map set to the zoomed extent.
I am using the ArcGIS Javascript API.
I have tried several JS/jquery/html methods to make this work, (window.load, ,$window.ready)..but nothing has worked for me.  I continue to get an error in my new page saying that 'map' is undefined, and I can see in Firebug that my function to set the new extent based on the extent of the original page runs before the page is fully loaded.
I hope I am explaining this problem clearly.  Does anyone have any ideas? Thanks,


texts: "Request payload size exceeds the limit" error when exporting image from Google Earth Engine

texts:
I have shape file with multi polygons and when i export though google earth engine this massage appear:

  "Request payload size exceeds the limit:41943304 bytes".

The problem was availability  of multi polygons in one shape file. Is there any procedure to solve the issue?

Here is the code:

Map.centerObject(AOI, 9)
var ndviParams ={min: 0, max: 1, palette: ['FFFFFF', 'CE7E45', 'FCD163', '66A000', '207401','056201', '004C00', '023B01', '012E01', '011301']};
var sentinal=ee.ImageCollection("COPERNICUS/S2")
var nir11 = image11.select('B8');
var red11 = image11.select('B4');
var ndvi11 = nir11.subtract(red11).divide(nir11.add(red11)).rename('NDVI').clip(A);
Map.addLayer(ndvi11, ndviParams, 'NDVI November');
Export.image.toDrive({image: ndvi11,description: 'NDVI_Nov',maxPixels: 3e12,scale:


texts: Why do I continue to get a datum conflict message even though layers match?

texts:
I have a document with multiple raster layers and a set of buffer zones for which I am trying to use the zonal statistics tool to summarize some of the raster information within the buffers. Most of the layers are in NAD 1983 and so I have projected my 30m DEM layer from WGS84 to NAD83 to match. However, I continue to get "001003: Datum conflict between input and output" warning message when I run the zonal statistics tool along with some other tools.
This is for the buffer layer
[1]: https://i.stack.imgur.com/pXlrx.jpg
This is for the DEM raster
[2]: https://i.stack.imgur.com/9TLUA.jpg
I am not sure if I am missing something obvious with the projection, and I have tried re-defining the projection multiple times.


texts: How can I contribute to openlayers examples list?

texts:
Who or where to can I post a request to add an example to the openlayers 3 examples list?

Particularly I would like to ask to add this simple custom control example because it is quite hard for new users to understand custom control class usage in the custom control example.


texts: Encryption of spatial data, GIS access and cloud synchronization

texts:
I'd like to encrypt a folder of many spatial datasets (rasters and shapes) because my office computer is shared with others. I've never used this kind of security but I think I should preserve my data from other guys.

So, my questions are:

If I encrypt a folder with password access, e.g. using VeraCrypt (suggestions?), how will GIS programs (QGIS, SAGA and R) handle internally opening files, geoprocessing and so on? Is this possible? Do I have to change permissions in the folder? What should I do to allow GIS programs to use that data and, on the other hand, keep safe from other people?
Is it possible to put that folder in cloud, to automatically upload the data and synchronize with another computer/backup?


texts: Discontinuous or Interrupted map projections?

texts:
Do any of the current open source GIS softwares support discontinous projection systems? E.g. sinusoidal, dymaxion, butterfly, etc? If not, are there any affordable commercial ones that do (affordable meaning < $150)?

Upate: apparently it was a bad idea to lead with "sinusoidal", as that projection also has a single point non-interrupted form. Please note, the key word here is discontinuous or interrupted. Thanks.


texts: How do I input this function to "Raster Calculator" in Model Builder , in ArcGIS 10.2?

texts:
The raster files that function refers are as follows:

Precipitation (o “Raster” in figure 1): Variable
CN_2S_7ma: Constant
CN_8S_7ma: Constant

The mathematical function that must be entered is as follows:

a) First condition: Condition for raster that have "No Value" or empty cells, they put a zero.
If "Precipitation"> 0, Response: Yes = going to Condition 2, No = 0

b) Second condition: If "Precipitation"> = CN_2S_7ma, Response:

Yes = 〖(Precipitation - CN_2S_7ma)〗^2 / (Precipitation + CN_8S_7ma)

No = 0

The idea is to place the function in a model builder, after using the "Iterate Rasters" to take all raster within the folder ("Precipitacion") and apply the function.
So far, I am writing like this both conditions in "Raster Calculator" and not working:

Con( “%Raster%” > 0,( Con( “%Raster%” >= "%CN_2S_7ma%",( ( Square( "%Raster%" - "%CN_2S_7ma%" ) ) / ( “%Raster%” + "%CN_8S_7ma%" ) ), 0 ) ), 0 )

Figure 1: print screen of my "not working" model Builder :(


texts: Creating boundary box in Python on randomly generated longitude / latitude points?

texts:
I'm trying to create a boundary box to separate a list of randomly generated longitude / latitude points. I would also like to look at the distribution of randomly generated points over time. I wrote a python script that randomly generates 5,000 longitude / latitude points.

How do I write the script and which libraries do I use?


texts:
I am running PostgreSQL 9.4 and PostGIS 2.3

For the background: I developped a DB on a regular machine used as host. I needed to use a view+rule system to allow the work to be done on spatial subsets of the datas. It worked fine, however when my company acquired a proper server, and after I migrated the DB, everything started to be extremely slow, a QGIS project is now taking around 2:30min to load, and the refresh timings are also extremely long.
I figure the problem comes from the performances allowed by the server, and since updatable materialized views are not possible, I would like to try and optimize my views query. (actually I am trying also something else: emulate a viewing system with tables and triggers, but it is a rather complex architecture).

I saw in this post PostGIS doesn't use spatial index with ST_Intersects that it is better to not put subqueries inside the st_intersects, which after testing obviously is conclusive. But I face another problem, my view is not only a subset of the raw datas, it also put in shape the foreign keys in order to be readable. So I need to use twice the same table: once to shape the foreign key, and another for the spatial subset. Here is an example of my original query:

SELECT object.id_obj,
    project.name,
    object.geom
   FROM object
     LEFT JOIN project ON project.id_project = object.id_project
  WHERE st_intersects((select geom from project where id_project = 1), obj.geom);

This query takes 123msec, with 61 results.

Here is what I tried:

SELECT object.id_obj,
        project.name,
        object.geom
       FROM object
         LEFT JOIN project ON project.id_project = object.id_project
      WHERE project.id_project = 1 and st_intersects(project.geom, obj.geom);

The query runs much faster, 12msec, but only get 46 results, because sometimes the objects have the same id_project as foreign key, but sometimes not, and the ones which don't are not included in the result.

I obviously see where the problem lies, but I cannot figure out how to solve it. Any idea?


texts: How to plot Geoserver data (point data) as image on the map instead of vector plot?

texts: Need the help to integrate it with angular 5 if possible.


texts:
I have two tables from Open Baltimore Dataset - Vital Signs (Housing and Community Development) and Housing Typology 2008. I am trying to find correlation between two geometries present in those datasets. Data information:

Vital Signs: 55 rows

Housing Typology: 710 rows

I want to find which row of Housing Typology belongs to Vital Signs. And both use different SRIDs as well:

SELECT Find_SRID('baltimore', 'vs14_housing', 'geom');
> 102685
SELECT Find_SRID('baltimore', 'housing_typology_2008', 'geom');
> 2248

And when I executed the following query:

select
a.gid as Vital_Sign_Index, b.gid as Housing_Typology_Index
from
    baltimore."housing_typology_2008" b, baltimore."vs14_housing" a
    where
    ST_Intersects(ST_Transform(a.geom, 2248), b.geom)
    order by b.gid

It returns something like this:

vital_sign_index  housing_typology_index
51  1
55  1
51  2
55  2
55  3
51  3
55  4
51  4
55  5
51  5
51  6

Put I want to do one on one mapping rather than many to one based on maximum intersection with that geometry. So in simple words, I want to group by based on maximum overlapping distance.

Currently it returns 1255 records due to overlapping regions while I want it to return 710 (Numbers of Housing Typology rows).


texts: What Spatial Reference System do I store Google Map's Lat/Lng in

texts:
I am new to PostGIS, Spatial Reference Systems and projections and want to store lat/lng coordinates retrieved from Google Maps and Openstreetmap and its services like the Geocoder/Nominatim. I'm working with data in America and the main calculations are for nearest neighbors. Results will be plotted back on a web map (eg: via Google maps API or leaflet)

Question: What spatial reference system is Google map using? I figured out so far that its the WGS 84 datum with the Mercator projection in the geographic coordinate system. Should I store the location data as it is, or transform to US National Atlas Equal Area EPSG:2163?

Since I will be calculating distances and finding nearest neighbors, if I guess correctly that doing nearest neighbor search requires transforming the entire table to EPSG:2163, then will the solution be to store the data in both Mercator and EPGS:2163?


texts: Piping data between PostGIS 'raster2pgsql' and PostgreSQL using Python subprocess?

texts:
My goal is to create a system that will take a GeoTIFF format raster image and programmatically import it into a Postgres table using PostGIS raster2pgsql. In order to get this done with minimal intermediate files, I am using the Python subprocess module to pipe the stream from the raster2pgsql command into the psql command directly.

Here is some relevant code:

def importRaster(rasterFileName, rasterId):
    rasterFileName = config.rasterImportFolder() + "/" + rasterFileName
    tableName = "raster" + str(rasterId)

    env = os.environ.copy()
    env["PGPASSWORD"] = config.dbPassword()

    gdal = subprocess.Popen(("raster2pgsql", "-s", config.rasterSourceProjection(),
            "-I", "../" + rasterFileName, "-F", "-t", "100x100", "-d", tableName),
            stdout=subprocess.PIPE, cwd=os.path.dirname(os.path.realpath(__file__)), env=env)

    psql = subprocess.Popen(("psql", "-U", config.dbUser(), config.dbName()),
            stdin=gdal.stdout, stdout=subprocess.PIPE, env=env)

    print gdal.stdout.read()
    gdal.stdout.close()

This solution has always been fairly unstable, but when running on my Windows machine it was bearable.  However, I recently migrated my test environment to a Debian based distribution on a Raspberry Pi 3.  I was expecting a decrease in speed/performance, but there was also a significant decrease in reliability, with the command failing nearly every time.

I started printing the GDAL process output, as you can see in the code, and I discovered a very sporadic collection of errors occurring from attempt to attempt.  Sometimes I will get something to the tune of "Premature end of WKB data", sometimes I will get "relation raster6 does not exist", and sometimes the piped data will completely miss the commit query at the end of the stream and no action will be taken at all.  It seems to me that the pipe from GDAL to PSQL is randomly skipping or corrupting the queries from the output SQL.  This would explain a missing or corrupted CREATE TABLE causing the "relation does not exist" error.  It would also explain corrupted WKB data strings and the missing commit statement.  If I run these commands manually in the console with the pipe, the process works completely fine.

I have considered running these processes separately and streaming the data from raster2pgsql into an SQL file as an intermediary.  However, these SQL files are quite large in my application (greater than half a gigabyte) and are not ideal for the types of machines this software will be running on.

Why does it seem that Python subprocess is dumping some of the SQL stream during the pipe?

I would also like to add that I am using a very similar process for piping from shp2pgsql to psql.  In this case, there is much less data due to smaller file sizes and the pipe works without error consistently.  It seems to be a volume issue, but I was under the impression that the subprocess module could handle large data volumes.


texts: Why are Export/Import Script, Edit, Debug greyed out from context menu of many tools in ArcToolbox?

texts:
Most of the options of Export/Import Script, Edit, Debug are greyed out on the ArcTools as shown below, did anyone have similar problem? Know why and what needs to be invoked?

I'm using ArcGIS 10.2 for Desktop.


texts: Adding MS Access Table as X,Y Points using QGIS?

texts:
Can you add an MS Access table as X,Y points in the same way you can Add X,Y data in ArcMap using:

File --> Add Data --> Add X,Y Data


texts: GRASS i.landsat.toar returns empty raster

texts:
I am trying to use the GRASS module i.landsat.toar.

Currently the files created are empty and it spits out:

So it seems that it is not picking up something?  The location is set to the extent of the landsat image.  Does anyone have an idea where I am going wrong?

Becky


texts: Change datasource of multiple layers from one db to other db in QGIS project?

texts:
I have a QGIS project with around 100 layers referring to a postgres database.
Since we're migrating to a new server, the database will be transfered too, but will be given a new IP. The connection between all layers in my QGIS project will be lost, due to a slight alteration in the source connection.
Is there a way to change the datasource of all layers, instead of the need to refer to this new db, layer by layer?
(Extra info: the way how the postgres database is structured isn't changing, the IP is the only thing being changed.)


texts: Is it possible to have an animated SVG as a vector layer feature?

texts:
Is it possible to have an animated svg as a vector layer feature?  I'm using D3 to make the animation. Here's my animated svg: https://jsfiddle.net/eflowbeach/opfkh0n9/1/

I thought I was onto something when I got this to work:

// Turn the svg into an xml string for use in OpenLayers
var svgstring = (new XMLSerializer()).serializeToString(d3.select(svg.node()).node());

where svg is the animated version, but it only captures a snapshot of the animation.

I get this into OpenLayers like this:

 new ol.style.Style({
                image: new ol.style.Icon({
                    src: 'data:image/svg+xml;charset=utf-8,' + encodeURIComponent(svgstring)
                })
            });

After I got this to work, it makes sense why it is not animated; the string will only capture the current state of the svg. So, is there another way to do this or do I have to look into putting this on a canvas?


texts: How to generate pgrouting source and target with OSM?

texts:
I'm wondering how to generate with the given lonlat in OSM PostGIS entities.

In pgrouting's website, they are telling that source and target are identifier of the vertex but not telling if this an osm_id for OSM database or any. My only given to generate is lonlat of the potential source/target. I'm using osm2po to create routable database.

Can anyone give me a clue, so that I can make my map more dynamic?


texts: Calculating Closest Distance Between All Combinations of Set of Polygons using ArcGIS for Desktop?

texts:
I'm looking to calculate all combinations of the closest distance 2 polygons in a set of polygons in ArcMap 10.2.2. I have an advanced license. I'm assuming it has something to do with the Proximity toolset (Near Tool maybe? I've never done a "Generate Near Table" using polygon edges, only centroids and points).

Ideally there would be a simple way to loop through all combinations of polygons without creating a Python script (Although I'll write one if I need to, I'd just rather not if there's an easy tool available in an ArcMap toolbox).

For example, if I have a set of 4 polygons, I would want to know the shortest distance between 1&2, 1&3, 1&4, 2&3, 2&4, and 3&4 (all possible combinations).


texts:
I downloaded a geotiff from the recently published 25m-DEM of the European Union (see: http://courses.neteler.org/eu-dem-new-digital-surface-model-at-25m/). One tile has about 5GB so I would like to crop it to a smaller extend. I'm trying to use gdalwarp but I can't figure out the correct commands ...

Here's my bounding box: N48 W15 S43 E5

This is the command I tried :

gdalwarp -of gtiff -t_srs EPSG:3857 -te 5 43 15 48 input.tif output.tif

This is the error message that gdalwarp returned:

  "Creating output file that is 0P x 0L. ERROR 1: Attempt to create 0x0
  dataset is illegal,sizes must be larger than zero."

Could anybody give me a hint?


texts: Updating 30 Shapefiles, adding attribute and populating based on data

texts:
I need to add an attribute to each shapefile attribute table in ArcMap, I was thinking to do it in SQL, but I'm not sure how to go about doing it since I'm very new to ArcMap.

Can someone explain to me how I would go about doing this?

The added attribute would need each value populated based on another attribute in the table.

Thanks.


texts: Double cascading value relation in QGIS attribute form

texts:
I run a PostgreSQL database using QGIS as a frontend. I have a table in which I store archaeological finds related to sites. The finds-table has no geometry and just stores date. To be able to query my data in different ways, I have three fields to describe the nature of these finds, which are however interdependent (material,group,type). E.g. material=Pottery, group=Vessel, type=Vipperow or material=Iron, group=Sword, type=A27. I have used the value relation widget to create cascading input forms, that should allow only object types from groups that can be e.g. iron to displayed, when I set material to iron and only types of swords, when I set group to sword.

These two relations look like that in the editor:

This worked fine as long as I only had one relation between material and group, now I entered the second one between group and type and this appears not to work. I created relation tables for both, which both look quite alike (see below) and use the same filters. However, the material-group relation works but the group-type relation does not. It shows certain possible input option but always sets the same in any case. It is also not possible to enter certain option, which are int the tables. The table I use for the relationships are.

Table with the relationship material to group = finds_groups

Table with the relationship group to type = finds_types

Is it possible to run double cascading value relations in QGIS, or is there some error in my settings?


texts:
Besides reading about bathymetry, I would like to work on some hands-on training that I can learn about bathymetry. Is there any good bathymetry out there that I can get some training with the datasets?
I want to learn how it works and is it a good for small locate area or wide area ? What would I like to know the pros vs cons about using it for the Marine areas ?


texts: Merging .asc files = "The command line is too long"

texts:
Having problems merging LiDAR .asc folder files to produce a single GeoTiff file.

I have followed various instructions e.g. the 2 methods of selecting input files for merging .asc files using

folder -> select all and
Choose input directory instead of files check box

Output folder & file created

When I click OK get - " the command line is too long " dialogue box error. Both input methods result in getting the same error message.

Any thoughts?


texts: Is it possible to make an installer for a published tool, to install the tool to someone else's ArcGIS server?

texts:
I have a geoprocessing tool that works on my instance of ArcGIS for server. Is there a way to easily install this on someone else's Server, without having them run the tool and publish it through ArcGIS desktop?

For instance, I know my code is in: C:\arcgisserver\directories\arcgissystem\arcgisinput\ServiceToolbox...

Obviously I can't just copy this directory structure over to another server, and expect it to recognize it and work. But that's basically what I want to do.


texts:
I have a bunch of linear data (trails) stored in PostGIS and visualized in QGIS for cartographic production.  Many trails have switchbacks, and are generally very "weavy", which makes for difficult and often inefficient labeling of trail names.  In many cases, the trails won't even label because it cannot find a suitable path, even when using the "curved" placement options:

Missing labels, too curvy:

I've tried simplifying (straightening) the geometry through a PostGIS view to create a sort of "label layer", but this seems inefficient (lots of layers, and doesn't scale well across zoom levels).

Is there a way to take advantage of the internal QGIS Label Expressions to simplify the geometry of the labels only (not the features themselves) to accomplish straighter, simplified labels and achieve a higher level of success and visual appeal for labels?

I was looking under Labels > Value > Expression Dialog > Geometry > Simplify but couldn't figure out a way to include this into the features "name" field in the label expression.


texts: WMS GetFeatureInfo request with buffer doesn't return any result (OpenLayers/MapServer)

texts:
I'm using OpenLayers 4 to build a GetFeatureInfo URL with a buffer parameter to a WMS layer in MapServer in order to find the closest feature of the layer given some specific coordinates.

The problem is that the request never returns any result at all. So I don't know if there is something wrong in my code or I'm missing something.

Here is my code explained:

var layer = new ol.layer.Image({
    source: new ol.source.ImageWMS({
        url: 'http://www.juntadeandalucia.es/medioambiente/mapwms/REDIAM_PPHH_2012?',
        params: {'LAYERS': 'MSPF_POLY'},
        serverType: 'mapserver',
        crossOrigin: 'anonymous'
    })
});

// I create the layer only for requesting purpose so I don't add it to the map.

var coordinates = [311618.514508171, 4040090.26715453];

var srs = this.map_.getProjection().code; //"EPSG:25830"

const getFeatureInfoParams = {
        'INFO_FORMAT': 'text/plain', //or application/vnd.ogc.gml
        'FEATURE_COUNT': 1,
        'SRS': srs,
        'Buffer': 10;
};

var viewResolution = this.map_.getMapImpl().getView().getResolution();
// 353.3212194629518

var urlToRequest = layer.getSource().getGetFeatureInfoUrl(coordinates, viewResolution, srs,
getFeatureInfoParams);

return urlToRequest;
// urlToRequest is ==> http://www.juntadeandalucia.es/medioambiente/mapwms/REDIAM_PPHH_2012?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&FORMAT=image%2Fpng&TRANSPARENT=true&QUERY_LAYERS=MSPF_POLY&LAYERS=MSPF_POLY&INFO_FORMAT=text%2Fplain&FEATURE_COUNT=1&SRS=EPSG%3A25830&Buffer=10&I=50&J=50&CRS=EPSG%3A25830&STYLES=&WIDTH=101&HEIGHT=101&BBOX=293775.7929252919%2C4022247.5455716513%2C329461.2360910501%2C4057932.988737409

I have seen that it always return I=50/J=50 and perhaps this could be the problem I don't know. It is that correct?

The pixel coordinates shouldn't be different in different bounding boxes values per request?


texts: pythonaddins interactive output registry to a shapefile recursively

texts:
I am using pythonaddin's onMouseDownMap to get flow paths interactively and add each line drawn to a shapefile. "line" shapefile is created each time using arcpy.RasterToPolyline_conversion whenever I click on map but sometimes not registered in AddasStreams.shp. Any suggestions where this is falling?

class FlowPaths(object):
    def __init__(self):
        self.enabled = False
        self.cursor = 3

    def onMouseDownMap(self, x, y, button, shift):
        # add line to view -- addOutputsToMap
        arcpy.env.addOutputsToMap = True
        point = arcpy.Point(x,y)
        ptGeometry = arcpy.PointGeometry(point)
        theLine = arcpy.sa.CostPath(ptGeometry, optfolder + '/dem', optfolder + '/flowdir', "BEST_SINGLE")
        arcpy.RasterToPolyline_conversion(theLine, "line", "NODATA", "", "NO_SIMPLIFY", "VALUE")
        searchcursor = arcpy.da.SearchCursor("line", ("OID@","SHAPE@"))
        firstrow = searchcursor.next()
        insertcursor = arcpy.da.InsertCursor(optfolder + '/AddasStreams.shp',("Id","SHAPE@"))
        insertcursor.insertRow(firstrow)

        del firstrow,insertcursor,searchcursor


texts:
Im looking foward to create a mobile leaflet application to collect data.

First i use the L.locate-control plugin to get my current location:

L.control.locate({
strings: {
    title: "current position"
},
locateOptions: {
           timeout: 10000}
}).addTo(map);

Then i initialize a button with the L.EasyButton plugin to store the current location with a marker and add it to my layerGroup (places):

var places = new L.layerGroup();

L.easyButton('fa fa-thumb-tack', function(){
map.addOneTimeEventListener('locationfound', function (event) {
var marker = L.marker(event.latlng, {icon: smallIcon2}).addTo(places);
});
}).addTo(map);

In the end im using :

 document.getElementById('export').onclick = function(e) {

        var data = places.toGeoJSON();

        var convertedData = 'text/json;charset=utf-8,' + encodeURIComponent(JSON.stringify(data));

        document.getElementById('export').setAttribute('href', 'data:' + convertedData);
        document.getElementById('export').setAttribute('download','data.geojson');
    }

to export the marker position as a geojson data.

Everything works well so far!

Now im looking for a way to edit the popupcontent of a marker by clicking the marker and adding a value/text with the keyboard. Similiar to the way it is possible with leaflet.draw.

How to achieve this?


texts: Indices similar to EVI for vegetation health analysis using UAV multispectral data

texts:
I know NDVI is the most commonly used index to start with vegetation health analysis. I am looking for other indices that can be calculated when the available UAV imagery has NIR, Red, Red Edge, Blue and Green bands. I found EVI, but since I  don't have satellite data, so that's out of the question as the C1 and C2 coefficients don't apply here. I can't even find MSI as I don't have the SWIR band. Similarly, I can't find CCI, because I don't have the TIR band. What other indices can I calculate using the bands I have, when the purpose is to assess whether the vegetation is stressed or not and how the photosynthetic activity is in a arid desert environment?
EVI = 2.5 * ((NIR – Red) / ((NIR) + (C1 * Red) – (C2 * Blue) + L))


texts: I'd like to create multiple layers for asset management using openlayers, could anybody help?


texts:
I understand that coordinate system mismatches between layers is a relatively common problem. Nevertheless, exploring other posts hasn't helped me understand my issue. In short, I have a DEM raster and a KML route, and they should overlap based on geography alone, but they do not.

DEM File

I have downloaded a Vermont statewide DEM raster (the top option, "ALL_Elevation...") from this link: http://vcgi.vermont.gov/opendata. The data is in the NAD83 coordinate system per the included metadata text file.

KML File

I also have a KML file of a hike in Vermont in the WGS84 coordinate system. I'd post it here, but unfortunately, I cannot find a method of uploading the file with this post! It was generated using a GPS watch, and it happily imports into Google Earth without an issue. (I can send along the route if it may help troubleshoot this issue. Please just let me know.)

Attempts at a Solution

I've tried to first add the DEM file, set the project coordinate reference system (CRS) to NAD83 to match the DEM, enable "on the fly" CRS transformation, and then add the KML file. This did not work.

The opposite order didn't work either (i.e., first add the KML file, set the CRS to WGS84, enable "on the fly" CRS transformation, then add the DEM).

I've also attempted to add just the KML file to a new project and "reproject" the file using Save As and saving the file with a new name and using the NAD83 CRS. Adding the new and "reprojected" layer didn't overlap the DEM either.

Lastly, I've read in another post that data using a WGS84 CRS will have negative easting values. I would tack on the negative sign to all easting values just to try this, but I'm uncertain how to access the raw KML data.

Software Info

I'm running OS X 10.11.4 with QGIS 2.4.0 Chugiak.

UPDATE

Although I still haven't figured out what is going wrongly here, I do think that the issue is related to QGIS being unable to identify all the pertinent CRS information from the DEM layer. I think that a solution may involve defining a custom CRS using the "right" Proj4 code that matches my DEM. Proj4 codes can be found here: http://www.spatialreference.org. I've tried a couple different custom CRS ideas, but nothing works yet...


texts:
I will be using mongodb database which supports data on an Euclidean plane (not spherical). I will be performing queries which return all points within a box.

For the box I need to specify the upper left and bottom right coordinates (http://docs.mongodb.org/manual/reference/operator/query/box/#op._S_box).

My question is based on an a lat/long coordinate as input, how can I create these two coordinates needed?

I want to take my coordiante and create a box of roughly ~1km (I understand since it is a box it will not be 1km in all directions, but I am not looking for perfect accuracy).

Is there an algorithm for this?


texts: Set automatically the name ID of many mosaics in Google Earth Engine

texts:
I'm using a code wrote here in Stack in the question Mosaicking a Image Collection by Date (day) in Google Earth Engine to generate mosaics of Sentinel-1 by date (day) with images that cover a region. But in the code, the new Image Collection don't have a ID for each mosaic, anything like the date of the images or the date of images interval that it uses. There is any way to do this?

My code is:

var start = ee.Date('2014-10-01');
var finish = ee.Date('2018-03-31');

var collection = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterDate(start, finish)
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
.filter(ee.Filter.eq('instrumentMode', 'IW'))
.filterMetadata('resolution_meters', 'equals', 10)
.filterBounds(poly);

print(collection,"Original collection");

// Difference in days between start and finish
var diff = finish.difference(start, 'day');
print(diff, "Diferenca em dias");

// Make a list of all dates
var range = ee.List.sequence(0, diff.subtract(1)).map(function(day){return start.advance(day,'day')});
print(range, "Lista de todas as datas");

// Funtion for iteraton over the range of dates
var day_mosaics = function(date, newlist) {
  // Cast
  date = ee.Date(date);
  newlist = ee.List(newlist);

  // Filter collection between date and the next day
  var filtered = collection.filterDate(date, date.advance(1,'day'));

  // get date as YEARMONTHDAY. For example, for January 8th 2010
  // would be: 20100108
  var date_formatted = ee.Number.parse(date.format('YYYYMMdd'))

  // make date band as an 32 bit unsigned integer and rename it as 'date'

  var dateband = ee.Image.constant(date_formatted).toUint32().rename('date');

  // Make the mosaic
  var image = ee.Image(filtered.mosaic());

  // Add the mosaic to a list only if the collection has images
  return ee.List(ee.Algorithms.If(filtered.size(), newlist.add(image), newlist));
};

// Iterate over the range to make a new list, and then cast the list to an imagecollection
var newcol = ee.ImageCollection(ee.List(range.iterate(day_mosaics, ee.List([]))));
print(newcol, "Collection with mosaic images");

Map.centerObject(poly);
Map.addLayer(newcol, {bands: 'VV', min: -20, max: -5});

Follow the link to my code: Mosaicking by date.


texts: Passing table from model to script in same model

texts:
In my case, i have a model, that produce 2 Tables (zonal statistis). I have a script that should work with that tables. How do I pass these tables to a Script? my script takes certain fields in these 2 tables and store them in new table (+ some operations). Model looks like this:

And the script looks like this:

import arcpy

result=arcpy.GetParameterAsText(0)

enviro=r"D:\Skola\Kart+geoinfo\Diplomka\Diplomka.gdb"
arcpy.env.workspace=enviro
new_T=arcpy.CreateTable_management(enviro,result)
arcpy.AddField_management(new_T,"ID","SHORT")
arcpy.AddField_management(new_T,"U_Vetrolamu","DOUBLE")
arcpy.AddField_management(new_T,"Od_Vetrolamu","DOUBLE")
arcpy.AddField_management(new_T,"rozdil","DOUBLE")
arcpy.AddField_management(new_T,"rozdil_procenta","DOUBLE")

with arcpy.da.SearchCursor(table,["MEAN","OBJECTID"]) as sCur:
 with arcpy.da.InsertCursor(new_T,["U_Vetrolamu","ID"]) as iCur:
      for row in sCur:
         iCur.insertRow(row)

data = {}
with arcpy.da.SearchCursor(table1, ["OBJECTID", "MEAN"]) as sc:
for row in sc:
    try:
     data[row[0]] = data[row[0]] + row[1]
    except KeyError:
        data[row[0]] = row[1]

with arcpy.da.UpdateCursor(result, ["ID", "Od_Vetrolamu"]) as uc:
     for row in uc:
        row[1] = data[row[0]]
        uc.updateRow(row)
arcpy.CalculateField_management(result, "rozdil", "[U_Vetrolamu] - [Od_Vetrolamu]", "VB", "")
arcpy.CalculateField_management(result, "rozdil", "([U_Vetrolamu] - [Od_Vetrolamu]/[U_Vetrolamu] - [Od_Vetrolamu])*100", "VB", "")

I need to use table and table1 from model in script...


texts: Iterative selection of features and export to shapefile using PyQGIS

texts:
I'm working with QGIS 3.10 A Coruña on MacOS 10.13 environment.
I have a vector polygon layer and I need to select features with a given attribute value and export them as a new shapefile. I'm a newbie in these matters and haven't get any further than
layer = qgis.utils.iface.mapCanvas().currentLayer()
request = QgsFeatureRequest(QgsExpression("CODNUT2 = 'ES24'"))
features = layer.getFeatures(request)

What should be the next steps to generate the shapefile with those selected features?
How to make it an iterative process to generate new shapefiles for every value in the CODNUT2 attribute?


texts: Don't understand error : Error in if (nchar(projargs) == 0) projargs <- as.character(NA) :

texts:
When I run two different scripts from examples (which seem to work where I took them from) I get an error in both cases. Something seems to be wrong in my general set up or something. But I have no idea where it's going wrong.
The error I get is the following:

    Error in if (nchar(projargs) == 0) projargs <- as.character(NA) :
    missing value where TRUE/FALSE needed

What does the error mean and how do I solve it?
As far as I understand something seems to be wrong in the projection part, but I believe this should work by just running the code, since it should be  working (they are both stand alone examples...)

The codes I run where I get the errors are:

library(sp)
set.seed(357)
pts <- data.frame(x = rnorm(100), y = rnorm(100), var1 = runif(100), var2 = sample(letters, 100, replace = TRUE))
coordinates(pts) <- ~ x + y ## This is where I receive the error
class(pts)
plot(pts)
axis(1); axis(2)

ply <- matrix(c(-1,-1, 1,-1, 1,1, -1,1, -1,-1), ncol = 2, byrow = TRUE)
ply <- SpatialPolygons(list(Polygons(list(Polygon(ply)), ID = 1)))
ply <- SpatialPolygonsDataFrame(Sr = ply, data = data.frame(polyvar = 357))
plot(ply, add = TRUE, border = "red")

(Taken from how to overlay a polygon over SpatialPointsDataFrame and preserving the SPDF data?)

and

r1 = cbind(c(180114, 180553, 181127, 181477, 181294, 181007, 180409,
180162, 180114), c(332349, 332057, 332342, 333250, 333558, 333676,
332618, 332413, 332349))
r2 = cbind(c(180042, 180545, 180553, 180314, 179955, 179142, 179437,
179524, 179979, 180042), c(332373, 332026, 331426, 330889, 330683,
331133, 331623, 332152, 332357, 332373))
r3 = cbind(c(179110, 179907, 180433, 180712, 180752, 180329, 179875,
179668, 179572, 179269, 178879, 178600, 178544, 179046, 179110),
c(331086, 330620, 330494, 330265, 330075, 330233, 330336, 330004,
329783, 329665, 329720, 329933, 330478, 331062, 331086))
r4 = cbind(c(180304, 180403,179632,179420,180304),
c(332791, 333204, 333635, 333058, 332791))

sr1=Polygons(list(Polygon(r1)),"r1")
sr2=Polygons(list(Polygon(r2)),"r2")
sr3=Polygons(list(Polygon(r3)),"r3")
sr4=Polygons(list(Polygon(r4)),"r4")
sr=SpatialPolygons(list(sr1,sr2,sr3,sr4)) # This is where I receive the error
srdf=SpatialPolygonsDataFrame(sr, data.frame(cbind(1:4,5:2), row.names=c("r1","r2","r3","r4")))

taken from the 'over' functions in the sp package help page.


texts:
I'm having a little trouble using the merge command in QGIS. All the files have a range of different attributes, with some attribute headers in common. Some of these have different data format types for the common column, so when I try to merge them I get an error message.

I'm happy to delete all columns other one I've batch inserted, but how can I batch delete all columns other than one? I know I can use the Drop Fields Tool (/Delete Columns) but I'd have to individually select each column I want to delete in each of the layers.


texts: Creating georeferenced random data in Python/R?

texts:
I'm writing a dissertation on different parameter settings in GWR . In order to perform tests I'd need the availability of geo-referenced data that show some sort of spatial correlation (i.e. values are not independent of their geographic position but vary according to space having zones with high values as well as zones with low values). Datasets should have geographical coordinates (i.e. latitude and longitude) and their density patterns should vary. As an example: data with one/two/three centroids or with uniform distribution or data concentrated more in a corner or on the edges of the map, etc..

A good starting point for me was reading  How can I conduct Geographically Weighted Principal Component Analysis using ArcGIS, Python and SPSS/R?

However I'm stuck in the process of giving these data a spatial autocorrelation.

Any idea on how to achieve this?

I seek code in either Python or R.


texts:
Is there a way to get a georeferenced profile graph from Stack Profile?

I have a DSM and use transects to get a profile graph from the DSM.
This gives an x,y value between 0 and 300m (the length of the transect).
I would like to get the profile graph in UTM.

Is this a possibility in Arc Pro?


texts: Buffer regions how do you get only a certain distance from the vector

texts: Currently have a few vector lines that i want to create buffers for e.g. 10m, 20m,30m etc. I have created buffers but its not what i was aiming for, I want to see how another vector layer is affected. For more specifics I want to have a buffer region that goes from 0-10, 10-20, 20-30 etc but currently i have buffers that start from the original vector so its ending up as 0-10, 0-20, 0-30. Is it possible to create buffers for certain distances away from the interested vector or not? If so, how can I do this?


texts: Html2canvas cannot render map in OpenLayers 3

texts:
I'm using OL3 and the last version of Html2Canvas (1.0.0-rc.5). When I try to make a PNG screenshot of the div that embed the map it seems that the map itself is blank. Also the feature, that are located on my map, are not drawn on the PNG produced by Html2Canvas.

html2canvas($('#mapDiv')[0]).then(function (canvas) {
        var img = canvas.toDataURL('image/png');
        var html = '<title>' + textHead + '</title>';
        html += '';
        html += '';
        var w = window.open();
        w.document.open();
        w.document.write(html);
        w.document.close();
    });


texts: Using raster calculator with GDAL in python?

texts:
i am using an hyperspectral image with 158 bands. I want to calculate raster. i am using the following python script

import gdal # Import GDAL library bindings
from osgeo.gdalnumeric import *
from osgeo.gdalconst import *
import pylab as plt
import numpy as np
import xlrd
# The file that we shall be using
# Needs to be on current directory
filename = ('C:/Users/KIFF/Desktop/These/data/Hyperion/10th_bandmathref')
outFile = ('C:/Users/KIFF/Desktop/These/data/Hyperion/Math')
XLS=('C:/Users/KIFF/Desktop/These/data/Coef/bcoef.xlsx')
wb = xlrd.open_workbook(XLS)
sheet = wb.sheet_by_index(0)
sheet.cell_value(0, 0)

g = gdal.Open(filename, GA_ReadOnly)

# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
    print ("Problem opening file %s!" % filename)
else:
    print ("File %s opened fine" % filename )

#band_array = g.ReadAsArray()
#print(band_array)
print ("[ RASTER BAND COUNT ]: ", g.RasterCount)

for band in range( g.RasterCount ):
    print (band)
    band += 1
    outFile = ('C:/Users/KIFF/Desktop/These/data/Results/Temp/Math_1_sur_value'+str(band)+'.tiff')
    #print ("[ GETTING BAND ]: ", band )
    srcband = g.GetRasterBand(band)
    if srcband is None:
        continue
    data1 = BandReadAsArray(srcband).astype(np.float)
    print(data1)
   # for i in range(3,sheet.nrows):
    b=sheet.cell_value(band+2,1)
    #print(b)
    dataOut = (1/data1)
    driver = gdal.GetDriverByName("ENVI")
    dsOut = driver.Create(outFile, g.RasterXSize, g.RasterYSize, 1)
    CopyDatasetInfo(g,dsOut)
    bandOut=dsOut.GetRasterBand(1)
    BandWriteArray(bandOut, dataOut)

for the print(data1) i got only some "1", but the real values are some floats

0
[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
1
[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]
2

  Pixel value   0,139200

Can you find the mistake?


texts: How to specify output table name in GeoPackage when using ogr2ogr

texts:
I want to convert a GeoJSON file to a Layer within a GeoPackage.

For this, I'm trying to use the following command:

ogr2ogr dist.gpkg mah.geojson -nlt PROMOTE_TO_MULTI

This sucessfully creates a new GeoPackage, with a layer whose name is mah.

I want to have control over the name of the output layer.

I have also tried adding a -lco IDENTIFIER=dist_2019 parameter, but this does set the name to what I want.

What is the correct way to specify the name of the output Table/Layer in the GeoPackage?


texts: Close gaps between some polygons, which are adjacent in reality

texts:
I have tens of thousands of polygons (farmer plots) which have been measured while on the fields. This has caused plenty of gaps between some polygons, which is not the reality. I would need to close these gaps, see the red marks (arrows and circles) in the screen shot attached.

However, I can't seem to find a good approach to only adjust the relevant polygons. I have mainly found solutions for when the polygons are supposed to cover an entire area, but this is not the case. I have tried the convex hull function but it is changing the shape of all the polygons, which I don't want.  Thus, I don't think that buffer and dissolve is an option. I also need to save the information kept in the attribute table. I only want the boundaries to meet, not overlap as I want to avoid double counting of some areas in my further analysis (yes, as also in the picture some plots overlap, which I need to go about, but that's another question).


texts:
I am writing a QGIS plugin in which one of the tasks is to fetch data (nodes, ways and relations) from open street map and display it on the plugin as a .shp file. Right now it returns the number of nodes and not the actual data.

My code snippet is as follows:

import overpy

class GetData():
    def __init__(self, p1, p2, console):
        print('getData Module Accessed')
        self.console = console

        api = overpy.Overpass()
        result = api.query('node({}, {}, {}, {}); out;'.format(p1.y(),p1.x(),p2.y(),p2.x()))
        nodes = len(result.nodes)

        consoleNodes = 'OSM Node:'+'<span style=\" color:#ff0000;\">'+'{}'.format(nodes)+'</span'
        self.console.append(consoleNodes)

self.console just shows the output.


texts: Using only part of inline variable substitution string in ModelBuilder?

texts:
I am using an inline variable substiution to alter a field name. I do this within an iteration model in model builder. The inline variable is the name of the feature class which I use as my input for each iteration. When I am altering the field name I would like to only use the first part of the inline variable. To explain this a bit better, here is the name structure of my feature classes:

SA_05_Walk_CBR
SA_10_Walk_CBR
SA_05_Bike_CBR
SA_10_Bike_CBR
SA_05_Driv_CBR
SA_10_Driv_CBR
SA_05_Walk_PJC
SA_10_Walk_PJC
SA_05_Bike_PJC
SA_10_Bike_PJC
SA_05_Driv_PJC
SA_10_Driv_PJC
etc.

(Note that all feature class names have the exact same amount of characters for each part of the names)

Within my iterator, these feature class names represent my inline variable substitution string (called "count_field" in my model, see picture).

When I use the iterator, I join all the SA_xx_xxxx with a point feature class where I use the last three letters of the feature class name as a wilcard (e.g. "star_CBR"). For my next iteration I do the exact same thing but then using "star_PJC" as my new wildcard. These iterations work fine without any problems and it is also the reason why my feature classes need to end with "*_xxx".

Next in my iteration, I alter a certain field name so that they can be differentiated within the table of the feature class. Currently I use this in-line variable subsitution to rename the field (e.g. "SJ_Count_%count_field%"). Also this process works fine. Here is a picture of the alter field window:

(Please ignore the explanation mark; this picture is pure for a demo)

The thing what I want to do, however,is that it only uses part of the inline variable to alter my field name. The last 4 letters of my inline variable (e.g. "_CBR") can be removed within my new field name. So a field name would be called, for example, "SJ_Count_SA_05_Bike"

I have tried in renaming my field name again but I haven't found a way to remove the last 4 characters through an iteration. Anyone any ideas how I could do this in an efficient way? Using python within the alter field tool does not seem to be an option (as seen from alter field tool window picture). Here is a picture of my model:


texts: Outputting colour associated with category in text format

texts:
I've written a plugin which randomly generates fill colours for differing categories. I'm wondering if there's a way to iterate through the category renderer and output the RGB colour as text.

I have so far an example of the creation of the categories:

def styleMap(self):
    categories = []
    for cat in range(0,10):
        symbol = QgsSymbol.defaultSymbol(self.activeLayer.geometryType())
        layer_style = {}
        layer_style['color'] = '%d, %d, %d' % (randrange(0, 256), randrange(0, 256), randrange(0, 256))
        layer_style['outline'] = '#000000'
        symbol_layer = QgsSimpleFillSymbolLayer.create(layer_style)

        # replace default symbol layer with the configured one
        if symbol_layer is not None:
                symbol.changeSymbolLayer(0, symbol_layer)

        # create renderer object
        if cat == 0:
                category = QgsRendererCategory("", symbol, "")
        else:
                category = QgsRendererCategory(str(cat+1), symbol, str(cat+1))
        # entry for the list of category items
        categories.append(category)
    renderer = QgsCategorizedSymbolRenderer(self.field_name, categories)

From here I'd like to create a new function which loops through all the features in the layer and returns the colour of the renderer in RGB to save to an attribute file.

So, if self.field_name == 1, return the colour associated with category 1 in text format.

How would one go about iterating through all the category data and accessing the colour information?


texts: lyr.datasetName works with "CURRENT" mxd but not with mxd path

texts:
I have some mxds in wich I want to set all the layers a layer description.
This works as far as I use the code with current mxd, but if I give a path it has no results, and I also don't get any error messages.

arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("C:/base/example.mxd")
for lyr in arcpy.mapping.ListLayers(mxd):
    lyr.description = "asd"


texts: Cursor returns object reference not set when reading Annotation Layer in ArcMap

texts:
I have an annotation Layer of Geometry Type Polygon. There are few hundred features in it and I need to shrink each of them to a smaller size (for example reducing them by 50%). I am trying to  use the following sample code but getting

object reference not set to an instance of an object

error at the line
pFeatureCursor = pGeoFeatureLayer.FeatureClass.Search(pSpatialfilter, false);

because of which I am not able to test and make sure if the logic written after this line is correct or not.
Is there sample code or a correction in given code I can use to achieve the desired result?
I am using C#, ArcGIS SDK for .NET version 10.2, .NET Framework 3.5 and Visual Studio 2010.
 protected override void OnClick()
        {
            try
            {
                IMxDocument pMxDocument = ArcMap.Application.Document as IMxDocument;
                IActiveView pActiveView = (IActiveView)pMxDocument.FocusMap;

            string AnnoLayerName = "";
            AnnoLayerName = "MainPipeAnno";

            IAnnotationLayer pAnnotationLayer = null;

            ISpatialFilter pSpatialfilter = new SpatialFilterClass();
            pSpatialfilter.WhereClause = "annotationClassID = 'Leader Line'";

            pAnnotationLayer = GetAnnotationLayer(AnnoLayerName);

            IFeatureCursor pFeatureCursor = null;
            IGeoFeatureLayer pGeoFeatureLayer = pAnnotationLayer as IGeoFeatureLayer;

            pFeatureCursor = pGeoFeatureLayer.FeatureClass.Search(pSpatialfilter, false);

            IFeature pFeature = null;
            IAnnotationFeature pAnnotationFeature = pFeature as IAnnotationFeature;

            while ((pFeature = pFeatureCursor.NextFeature()) != null)
            {
                IGeometry pFeatureGeom = pFeature.ShapeCopy;

                // Create a point to act as origin at the centre of this polygon or collection of polygons using IArea interface
                IArea pFeatureArea = (IArea)pFeatureGeom;
                IPoint origin = pFeatureArea.Centroid;

                //Cast the feature geometry into the correct Interface
                ITransform2D transformScalePoint = pFeatureGeom as ITransform2D;

                //Perform Scale
                transformScalePoint.Scale(origin, .5, .5);

                // now the object stored in pFeatureGeom is reduced to 50% of size scaled around the centre of the feature
                pFeature.Shape = pFeatureGeom; // update the features' geometry

                pFeature.Store();
                pActiveView.Refresh();
            }

        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.Message);
        }
    }

    //This Function fetches the Annotationlayer from ArcMap depending upon Layer Name passed
    public static IAnnotationLayer GetAnnotationLayer(string strAnnotationLayerName)
    {
        IMxDocument pMxDocument = ArcMap.Application.Document as IMxDocument;
        IMap pMap = pMxDocument.FocusMap;

        IAnnotationLayer retAnnotationLayer = null;

        for (int i = 0; i < pMap.LayerCount; i++)
        {
            ILayer tempLayer = pMap.get_Layer(i);

            if (tempLayer.Name == strAnnotationLayerName)
            {
                retAnnotationLayer = (IAnnotationLayer)tempLayer;
                return retAnnotationLayer;
            }
        }

        return retAnnotationLayer;
    }


texts:
Might it be that the SQL API version available on GitHub is still v1.x.x?

See also https://github.com/CartoDB/CartoDB-SQL-API/blob/master/package.json

I mean, on cartodb.com they are using v2 already...


texts:
Possible Duplicate:
  Shapefiles to Text

I have several shapefiles from which I need to extract all of the points, with all of their attributes, including position coordinates. I only have access to opensource solutions, but can use whatever would be quickest/easiest, and would not be adverse to writing a little code if necessary. What would be the best way to go about this? Ideally I'd like the output to be a CSV.


texts: OS Master Map Styling in QGIS using QML/SLD

texts:
I am  trying to apply stylesheet that I downloaded in qml format from GitHub of Ordnance Survey https://github.com/OrdnanceSurvey/OS-VectorMap-District-stylesheets/tree/master/GML%20stylesheets

Master Map looks fine in the default styling but when I apply the style that I downloaded makes the Map disappear. What could be the reason to make it disappear into QGIS.

Same happens with SLD stylesheet

Below is the screenshot of the attributes of the QML stylesheet


texts: Setting Field in ArcMap Symbology Properties using ArcObjects?

texts:
I am programmatically rendering a feature layer in ArcMap but I can't set the Field property in the Symbology properties form. This results in the layer having the correct render/breaks/legend settings but the on screen render is wrong. When I open the Symbology tab the Field isn't set. When I do select the field the classification is reset. The focus field is actually from a linked table. Any suggestions? Below is the code I have for rendering the layer.

I'm using ArcGIS 10.2.2; Windows 7;VS2010 Express. Feature layer is in a file geodatabase along with a linked table.

I've asked this question on the ESRI forums with no solution yet so hopefully someone here can help.

    public void ClassifyLayer(ref IFeatureLayer pFL, ITable pTable) {

    IGeoFeatureLayer pGFLyr = (IGeoFeatureLayer)pFL;
    object dataValues;
    object dataFrequencies;
    ISimpleLineSymbol pLineFillSymbol;
    IClassBreaksUIProperties pClassBreaksUIProperties;
    IFeatureRenderer ifr = pGFLyr.Renderer;
    IClassBreaksRenderer pClassBreaksRenderer = (IClassBreaksRenderer)ifr;
    IDataSampling ds = (IDataSampling)pClassBreaksRenderer;
    ITableHistogram pTH = (ITableHistogram)new TableHistogram();

    pTH.Field = "IDX";
    pTH.Table = pTable;
    ds.MaxSampleSize =  pTable.RowCount(null) + 1;;
    pTH.Sampling = ds;
    IHistogram pHist = (IHistogram)pTH;
    pHist.GetHistogram(out dataValues, out dataFrequencies);

    IClassify pClassify  = new QuantileClass();
    double[] classes;

    pClassify.SetHistogramData(dataValues, dataFrequencies);
    pClassify.Classify(10);
    classes = (double[])pClassify.ClassBreaks;
    pClassBreaksRenderer.BreakCount = classes.Length - 1;
    pClassBreaksRenderer.Field = pTH.Field;
    pClassBreaksRenderer.MinimumBreak = classes[0];

    //Spectrum-Full Bright
    IMultiPartColorRamp pMultiColourRamp = GetSpectrumRamp(10);
    bool rampcreated = false;
    IEnumColors pColors = pMultiColourRamp.Colors;

    for (int i = 0; i < pClassBreaksRenderer.BreakCount; i++)
    {
    pLineFillSymbol = new SimpleLineSymbol();
    pLineFillSymbol.Color = pColors.Next();
    pClassBreaksRenderer.Symbol[i] = (ISymbol)pLineFillSymbol;
    pClassBreaksRenderer.set_Break(i,classes[i + 1]);
    pClassBreaksRenderer.Label[i] = ((Single)(classes[i])).ToString() + " - " + ((Single)(classes[i + 1])).ToString();
    }

    pClassBreaksUIProperties = (IClassBreaksUIProperties)pClassBreaksRenderer;
    pClassBreaksUIProperties.Method = pClassify.ClassID;
    pClassBreaksUIProperties.ColorRamp = "Spectrum-Full Bright";

    for (int i = 1; i < pClassBreaksRenderer.BreakCount; i++)
    {
    pClassBreaksUIProperties.set_LowBreak(i, pClassBreaksRenderer.Break[i - 1]);
    }

    pClassBreaksUIProperties.ShowClassGaps = true;

    pGFLyr.Renderer = (IFeatureRenderer)pClassBreaksRenderer;
    IRendererPropertyPage pRendererPropPage = (IRendererPropertyPage)new GraduatedColorPropertyPage();
    pGFLyr.RendererPropertyPageClassID = pRendererPropPage.ClassID;

    pMxDoc.ActivatedView.PartialRefresh(esriViewDrawPhase.esriViewGeography, pFL, null);
    pMxDoc.UpdateContents();

}