texts: Why extract returns NAs on the boundary when polygons exactly overlap with the rasters in R?

texts:
I'm new to GIS and English is not my first language, so apologies for anything inappropriate.

I'm trying to use extract() to get mean values of rasters inside polygons. One of the R codes I used are

county_NDVI<-extract(NDVI,province,fun = mean,df = TRUE,weights = TRUE,normalizeWeights = TRUE)

The CRS's of 'NDVI' and 'province' were different, but I've changed the CRS of 'province' with spTransform so they are the same now. And I have plotted them and checked that they are exactly overlapped. Also the raster data are continuous, so there are not supposed to be NAs.

However, out of 2896 polygons, I got 365 NAs. And I identified the location of all the NA polygons and found that they are all on the bounday of the polygons. So I suspect that the problem has something to do with the boundary or the units of the two objects (NDVI and province) are slightly different.

Any suggestions on what else I can check or what I can do to solve the problem?


texts: Exporting groups of PDF based on Orientation

texts:
I am trying to create a map book with two different orientations. The part I am currently stuck on is for the loop to go through the table and export groups of the same orientation. an example of the table is like:

_OR|Page_Num
Portrait|1
Portrait|2
Portrait|3
Landscape|4
Landscape|5
Landscape|6

and would like to have it exported as "Portrait 1-3.pdf"(I have another mxd for landscape so I would just make a quick edit so it works for landscape instead of portrait) so that I could combine the smaller groups in sequential order using bluebeam. I currently am using code I found on here that selects all portrait or landscape and exports them, but would be ideal if I could break it up.

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
ddp = mxd.dataDrivenPages
indexLayer = ddp.indexLayer
output_folder = r"C:\Users\JAlmeida\Downloads\Pdf"
select = set(row[0] for row in arcpy.da.SearchCursor(indexLayer, '_OR'))

for _OR in select:
    where_clause = '"_OR" = ' + "'%s'" %_OR
    arcpy.SelectLayerByAttribute_management (indexLayer, "NEW_SELECTION", where_clause)
        ddp.exportToPDF(str(_OR) + ".pdf", "SELECTED" )
del mxd


texts: Stop automatic ANSI_1252 to UTF-8 encoding when adding a DXF layer?

texts: I have a DXF file encoded in ANSI_1252, with object attributes containing umlauts. When I open it in QGIS, the umlauts don't display properly. This apparently happens because QGIS automatically sets the character encoding to UTF-8. In Layer Properties / General the Data source encoding is set to UTF-8 and greyed out, so I can't change it there either. How can I keep the ANSI_1252 encoding?


texts: Which interpolation method for noise mapping?

texts:
I want to create a noise map with measured point data. But I want to have a planar representation, so I need an interpolation method.

Which one is the 'best' for measured noise data?

I prefer 'IDW' interpolation, but I am not sure about that. Any ideas?


texts: Looking for collaborative webmapping projects

texts:
I´m looking for open source collaborative web mapping projects which allow vector editing, wms display and other features like comments, photo-uploads etc.

I´ve only found Ushahidi and Sahana so far, but I´m sure there must be others out there.

If anyone has heard of anything, please let us know.


texts: What is the faster way to run an OSM (xml) data to Esri geodatabase data for a big area with FME?

texts:
I am working on OSM highway(polyline) data with respect to the keys Max_speed and surface for the whole Africa. In this case I took at first the  Geometry Filter with line and I used 2 list_searchers for max_speed and surface.I also use 2 attribute_creator for output. This way is working to get a Geodatabase file which i can use on arcmap. But it takes a long time to read the data.

Can anyone help me to give the faster way to run the XML data? OR
the faster way to use PBF format in this purpose?


texts: QGIS - plot points on a map proportionally to column values

texts:
I have a .csv file like as follow:

latitude | longitude | norm_imp
-22.13   |  -43.18   |  0.017
-22.19   | -43.29    |  0.038
-22.33   |  -43.08   |  0.117
-22.19   | -43.29    |  0.2038

I plot the latitude and longitude points on a layer, but I need these points to be proportional to those in the norm_imp column. Is it possible to do this in QGIS? If so, how?


texts: SharpMap Coordinate transformation Problem

texts:
I want to transform my layers from SRID 31463 (Gauß Krüger) to Webmercator. Since SharpMap does only know webmercator and WGS84 i tried to create the Coordinate System 31463 by WKT and applied it on the layer. But i have a small offset of a few meters.

        ProjNet.CoordinateSystems.CoordinateSystemFactory csFact = new ProjNet.CoordinateSystems.CoordinateSystemFactory();

        string wkt_gk = "PROJCS[\"DHDN / 3-degree Gauss zone 3 (deprecated)\",GEOGCS[\"DHDN\",DATUM[\"Deutsches_Hauptdreiecksnetz\",SPHEROID[\"Bessel 1841\",6377397.155,299.1528128,AUTHORITY[\"EPSG\",\"7004\"]],AUTHORITY[\"EPSG\",\"6314\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4314\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",9],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",3500000],PARAMETER[\"false_northing\",0],AUTHORITY[\"EPSG\",\"31463\"],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH]]";

        GeoAPI.CoordinateSystems.ICoordinateSystem gauss_krueger_3 = csFact.CreateFromWkt(wkt_gk);

        GeoAPI.CoordinateSystems.ICoordinateSystem webmercator = ProjNet.CoordinateSystems.ProjectedCoordinateSystem.WebMercator;

        ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory ctFact = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();

        layer.CoordinateTransformation = ctFact.CreateFromCoordinateSystems(gauss_krueger_3, webmercator);
        layer.ReverseCoordinateTransformation = ctFact.CreateFromCoordinateSystems(webmercator, gauss_krueger_3);


texts: Moveable ArcGIS Identify Popup Window (Dojo)

texts:
I'm trying to get a popup window that identifies a location to be moveable. My code is similar to the API example:

var popupOptions = {
    fillSymbol: new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID,
                new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_SOLID,
                new dojo.Color([255, 0, 0]), 2),
                new dojo.Color([255, 255, 0, 0.25]))
}

var popup = new esri.dijit.Popup(popupOptions, dojo.create("div"));

map = new esri.Map("map", {
    infoWindow:popup,
    extent: initExtent
});

I can't seem to find anywhere that will allow me to make that popup infowindow moveable (was thinking of trying to implement dojo.dnd.moveable).

I've tried, unsucessfully:

function makeMoveable(node){
    var dnd = new dojo.dnd.Moveable(dojo.byId(node));
}
makeMoveable(popup);

Anyone have a suggestion?


texts: What is the "clip" feature for color ramps?

texts: In the singleband pseudocolor band rendering, there is a Boolean checkbox "clip". What does this do?  I don't notice a difference in my maps and haven't found documentation for it.


texts:
I'm creating a new WFS-T feature but I'm getting the below exception when I try to send the feature:

  No such feature type http://IP:PORT/geoserver/CARBONITA:bairros

link to gist:
 https://gist.github.com/oiacrasec/e82a0fe6e30cfd5f8770fbfd43b6a49e

gist contains:
 - Geoserver capabilities
 - XML sent
 - XML received
 - Project files


texts: Error message trying to merge Shapefiles in Python

texts:
I try to merge all Shapefiles from one folder together into one, using this:

shapeList = arcpy.ListFeatureClasses()
print shapeList
for fc in shapeList:
    if fc == "*domain*":
         arcpy.Merge_management(shapeList, "MergeList")

But I always get this error message:

Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 3987, in Merge
    raise e
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Datasets:
-- HERE THE LIST OF FEATURECLASSES --
does not exist or is not supported
Failed to execute (Merge).

I am sure there are only feature classes in the folder, when I try to use a Merge function I get the Error:

ERROR 000338: Inputs must be either all Feature Classes, Tables or Rasters; not mixed.
Failed to execute (Append).
Failed at Thu Jan 29 16:22:28 2015 (Elapsed Time: 0.00 seconds)

Does anybody know what's going wrong?


texts: Why does ModelBuilder lose data sources when copy/pasting or renaming toolbox?

texts: I have a ModelBuilder with a parent model that contains several sub-models.  The parent and sub-models are contained w/in the same toolbox.  Everytime I copy&paste or rename the toolbox, the parent models shows a red x and when i open it, I can see that the submodels w/in the parent also have red x's.  The issue seem to be that the sub-models loose their input sde data sources and so I have to re-establish them and then re-add them to the parent model.  I would like to find out what causes the models to fail like this when copying or renaming.  I'm using Desktop 10.2.  The screenshot is of the parent model.  The top row shows the failed sub-models and the bottom row is me in the middle of editing and re-adding the submodels to the canvas.


texts: Obtaining pixels between 2 raster layers (border)?

texts:
I have a rasterStack obtained from a stack of 3 raster layers by:

    image = stack(A, B, C)

I would like to calculate the distance between the pixels of layer A and pixels of layer C. After obtaining SpatialPoints of these 2 layers I used:

    A_pts = as(red,"SpatialPoints")[red[]==1]
    C_pts = as(blue,"SpatialPoints")[blue[]==1]
    Distance_AC <- gDistance(A_pts, C_pts, byid=TRUE)

But as the number of pixels are huge the calculating process cannot finish.

Therefore I would like to calulate the distance between the layer C and only the pixels of the layer A that are in contact with layer C = the distance between pixels of layer A and the pixels of the border between A/C.

How could get these pixels between layer A and C as illustrated in the  following image ?


texts: Selected features tool don't work in QGIS 2.12

texts:
I work with this polygon layer and features 6 and 7 are identical:

I don't understand why the selected feature tool work when i select feature number 7

and don't work when i select feature number 6?

Note: i copy/ paste feature 7 from feature 6
I red How to select features of visible layers in QGIS and that didn't solved my problem.


texts: How to Align Vector Layer on QGIS cloud Imageries for Printing

texts:
I noticed a weird problem every time I save maps with either Bing or Google Imageries from QGIS Cloud. The vector features all align perfectly well on the map canvas and even on print composer. But after saving the map as either .jpg or .pdf, the vector layer on the saved map file would have shifted by a huge margin.

Please see picture attached.

Does anyone have a solution?


texts:
I have found a very interesting paper here related to this topic.

I just want to have views/ideas from other professional about how we can protect GIS data using Cryptography or Digital Watermarking?

Any suggestions, especially for Administrative layers of a country where source of layers is same for every company so output is also same more or less.

I am interested in some methods like one suggested here


texts:
I have been working for some years with Google Maps and am now dipping my toes into the ArcGIS world. I am looking for some equivalent to Google's static maps api or Mapbox's  static image API which can capture an image from a particular built-in basemap (or perhaps custom map) at specified coordinates with specified dimensions and a small handful of additional overlays (single marker, polygon etc). This API would then return an image that I could embed in my web application via URL. And would have some mechanism for dealing with access tokens etc.
I am fairly new to the whole world of ArcGIS Online / ArcGIs Platform / ArcGIS Runtime / World Geocoding etc., so partly I'm just not sure where to look but so far in my searching I have not turned up anything.
Again, this is not for offline use or a "download" per se, just looking for a way to generate some simple images that I can embed in a website.
Does Esri / ArcGIS have such a service? Better yet, is it accessible via their Java APIs?


texts: ADD GIS functionality to existing ASP.Net tool/application

texts: I have an ASP.Net application that receives IOT data from field devices associated with clients. Currently we use a MySQL database for storing everything and much of the metadata associated with the devices is entered manually via forms. I've tried QGIS and pyGIS but am unclear as to how to incorporate into an existing application and how to satisfy expectations relating to incoming changes via GIS files. For example - If a client updates a pipe diameter on their shapefile and imports their updated file into the application, is there a tool that can bind to the database to look at the changes and incorporate those changes?


texts:
I was using Wroclaw 1.7/GDAL 1.9(earlier release)/Grass 6.4/GSL 1.15-2 and started getting the following message when saving any item to esri format.

"Export to vector file failed.
Error: creation of layer failed (OGR error:Failed to open Shape DBF file `/Users/jes....dbf)"
And when I look for the file/sibling files, the dbf is not there.

I just updated everything QGIS 1.8/GDAL 1.9 (December, 3 2012)/Grass 6.4.2. Still getting the same message with any feature that I save.

What can I do to fix this asap?


texts: arcpy.importDocument of ArcGIS Pro not working on specific MXD and not throwing error

texts:
I have a script that imports every MXD in a folder to an empty ArcGIS Pro project, then use saveACopy to save to a new project file. The problem is: The script just stops working every several MXD. Nothing happens, no error thrown.
For debugging, I got 2 print statements, one immediately before, one immediately after importDocument statement:
print(mxd_path)

aprx.importDocument(mxd_path)

print('imported')

Whenever the script reached a troublesome MXD, I got its path printed (the path looks totally normal to me), then nothing happened (have left it just like that for an hour). I've tried to import just that troublesome MXD with the script and got the same result: just the path printed, then nothing happened.
When I open the ArcGIS Pro app and import the manual way, i.e. right-click on the mxd file in Catalog pane and select Import and Open, the same troublesome MXD is imported and opened in seconds.
For MXD that can be imported via the script, it's done in less than one second.
MXD files are located in local hard drive. MXD file names contain A-Z character and dot. I've tried removing the dot from the file name of the troublesome MXD and that did not help.

The good news: When I opened the troublesome MXD in ArcMap application, saved it without doing anything, then re-run the script, that MXD was imported in less than one second. That means opening the MXD in ArcMap, then just save it will solve the issue.
The bad news: the way ArcPy opens and saves map Document is different from the way ArcMap application does. So using a script to open and save (or saveACopy) does not solve the problem.


texts: How to stay productive and proactive at work?

texts:
I know this is technically seeking advice, but, to me, advice on the subject is considered an answer to the question. Therefore, I wouldn't consider this to be outside the bounds of this StackExchange. I am looking for an answer that would include advice of what I should strive for and why it is important.

Context
I am new to GIS. I started a job doing part-time work for a place that had just launched a GIS side project. That side project is rapidly evolving into something bigger and I am the sole person in charge of making it happen. I didn't go to school for this. I am self-taught and learning as I go to provide what is requested.

Question
How do I continue to be productive and proactive at work?

Since I have no prior experience, I do not have many fresh ideas that I can present to help grow and evolve this project. Where can I seek inspiration in order to propose new projects and train myself? The last thing I want is to be working here and have everyone ask "...What does he even do?". My job includes both gathering data in the field and presenting it via printed maps. I use the ESRI ArcGIS suite.


texts:
Cell Statistics in ArcGIS will calculate the mean (or statistic of choice) among multiple rasters on a cell-by-call basis.
Can QGIS do this and, if so, how?


texts: How to export image on google earth engine for large area?

texts: The area that i have been work have about 159695,938 km² and I am mapping the deforestion on this with NDVI. I had to use 2 years to remove all the clouds (eg. 2015-01-31 , 2016-12-31) using a reducer.mean and cfmask band (surface reflatance Landsat 8 image), and when i export.to.drive i got the images separated. Is it possible to join this mean NDVI image in a single image to export? I will use this image in a geoprocessing software.


texts: How can I visualize a series of points on a world map?

texts:
I'd like to run a big list of IP addresses through geolocation and then display them on a map, most preferably an OpenStreetMap one. While geolocation will probably be easy using GeoIP, I'm looking for a free offline project that would allow me to visualise the GPS coordinates. How could I do that?

The platform I am interested in is Linux - it would be perfect if I could perform this task as a batch job so I could integrate it into my project: the input would be a list of GPS points, the output - an image showing hte map along with those points.


texts: display of raster in QGIS - cells are not rotated

texts:
when viewing rotated rasters in QGIS - the overall raster is rotated (as designed), though each grid cell is not. (the nodata value doesn't appear to pick up as the default in QGIS either - -2.9999e+30 instead of -3e+30, but that might be a separate question).

.

and a closeup

is this common for QGIS? viewing the same in 'the other GIS GUI' appears to rotate each cell.

i can add all the code - but the basic part is

def writeRast(gridbound, outTiff, grid, cellsize=1000, wkt=''):
    ''' lots of other stuff'''

    geotransform=(ulx, xres, xtheta, uly, ytheta, -yres)
    drv = gdal.GetDriverByName('GTiff')
    ds = drv.Create(outTiff, ncol, nrow, 1 ,gdal.GDT_Float32)  # Open the file
    band = ds.GetRasterBand(1)
    band.SetNoDataValue(-3e30)
    ds.SetGeoTransform(geotransform)  # Specify its coordinates
    if wkt != '':
        targproj = osr.SpatialReference()
        targproj.ImportFromWkt(wkt)
    band.WriteArray(grid)

is there some setting in QGIS that would display each cell rotated? or something i'm missing in writing the raster with GDAL?

the geostransform params (for at least of the rasters) look something like:

(-247544.61184672063, 900.39636198453309, 1333.8944867046337, 694821.1191712627, 1333.8944867046337, -900.39636198453309)

tried with gdal 1.11.0 and 2.0.
on mac OSX, QGIS v 2.8.3.


texts: Generating random points in selected areas?

texts:
I think my question is answered in different posts but I can't really make sense of them.

Background: We are designing a survey of a sample population of farmers to identify the seeds they use. I have a vector layer indicating the land use (irrigated or rainfed). This layer is composed by 16000 polygons spread across the country.

What I need: I want to randomly generate points within this polygons where to survey. More specifically, I want to generate 66 points in every province (provinces are indicated by another vector layer). In other words I want to generate 66 points under two conditions: 1) within province boundaries, 2) within the land use polygons.

I already tried QGIS plugin "Random points" but it seems to crash the program...

How should I proceed?


texts:
I want to use the tool intersects in the field calculator. If the condition is true, I want to set a result field to a specific number. The field calculator doesn't accept my expression but I can not find the reason for the error. Can you please help me.

CASE
WHEN intersects('layer1', 'layer2') = true THEN 1
WHEN intersects('layer1', 'layer3') = true THEN 2
END

the error says:

  syntax error, unexpected NUMBER_INT, expecting WHEN or ELSE or END
  syntax error, unexpected END, expecting $end

I work with QGIS 3.4 (Madeira)


texts: How to rotate QGIS composer map frame to North when using Atlas?

texts:
In QGIS, I'd like to use Atlas in the Print composer to generate a series of maps.  However, I'd like each map to be printed with "north up".  For example, using the classic projection of the United States (shown below), Washington and Maine end up being "tilted" so that north is not "up" on the page (see example).

Is there a setting in QGIS that will sync/rotate the map canvas to "true north" based on the current "zoom"?  I know I can create a new field with a "rotation" value and apply that to the properties of the data frame, but I've got too many features to calculate and populate manually.


texts: How could I classify every county in this US as flat, hilly, or moderate?

texts: I'm looking for a VERY rough label for every county, so low-res data is fine here. Is there data on existing elevation change per county? Or should I get an elevation layer of the US and do a zonal summary of intersecting features? Where could I find a rough elevation raster of the US?


texts: Using ArcGIS functionality in Python without running ArcMap (e.g. separate editor/environment)?

texts:
I'm trying to take the next step in programming ArcGIS with Python and see if its actually possible to run my script/code without having to load ArcMap and paste them into the Python window. I imagine I would do this view some sort of IDE specific to Python that loads something from ArcGIS.

How would I got about doing this?

Update: Summary of the solution

I went the Eclipse + PyDev route. Here are my notes on installing everything and getting it working:

Download Eclipse

This is slightly confusing given all the different flavors that are available (Java, Java EE, C/C++, etc.). I went with C/C++ since that's what familiar to me. It doesn't really matter since PyDev gets installed after this step. Pick what you want.

Install PyDev

PyDev will be installed via Eclipse. The instructions linked worked perfectly.

Configure PyDev

The above question gives a screenshot of a Windows configuration, while this link from an ESRI blog has more info for Unix systems although it is from 2008.

Ran a test script

Fired up a script I had been working on and it worked. There seems to be some differences in executing code in this way in how to access shapefiles and the workspace, etc, more explicitness is required. More to learn here.


texts: How to stop TerriaJS from automatically animating time-dependent data

texts:
I love TerriaJS, but one annoying feature is that with the default deployment, any time-dependent layer (like a time-dependent WMS) will start animating immediately.

Since some of our WMS layers take a few seconds to load each time step and have thousands of time steps, this is really more than annoying.

Is there a way to turn of this default behavior and allow the user to animate only when they click the play button on the time slider (at the bottom of the screen grab here):


texts: Counting points within range using CartoDB?

texts:
I have a dataset of users location and want to assign per user how many other users they have within a given range, and updating this value into a table.The dataset name is "users_20160410" and the column to be updated is "countdistance"

I have manage to get the number of total users within 10km for a given point, which is a unique value:

UPDATE users_20160410 SET countdistance = (select count(*) from users_20160410 where st_distance(the_geom,st_GeomFromText('POINT(-3.703550 40.417192)', 4326), true) < 10000)

What I want is to know for every single user how many other user he has within a range of 10km (this means get a value per user, not 1 unique as per my previous query). Also as the distance is not a fixed given point, the st_Geom_From_Text could not be used as it should take the exact location per user.


texts: Plotting shapely polygon with holes does not plot all holes


texts: Replace occasional Space with Underscore using Python parser of ArcGIS Field Calculator?

texts:
I am creating a model which will give me two string values, one which reads exactly what is in a field and another which automatically replaces any Space with an "_". I have used Calculate Value tool, and entered in the following:

Expression:

ReplaceName(%Value%)

Code Block:

def ReplaceName(name):
 return name.replace(" ", "_")

The result is an error:

ERROR 000539: Error running expression: ReplaceName(NoSpace1)
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
NameError: name 'NoSpace1' is not defined

I suspect because not all my Values have spaces

How can I make this work?


texts: ArcMap 10.2 crashes when running a Model

texts:
We have a feature class "A" with features that are modified by our users on a daily basis. One attribute is the "Date of Modification" that is being written for each edit done on a feature. We have another feature class "B" that contains the same features (with some different attributes) in which we want to "copy" all features that have been modified in feature class A and, at the same time, delete the same features in feature class B so that only the new modified ones exist.

We have created a model for this. It starts with the Select Layer By Attribute tool to select all features that have "Date of Modification" = CURRENT_DATE from feature class "A". Another Select Layer By Attribute tool selects the same features from feature class "B" (based on a common field, unique for each feature). Then the Delete Features tool is used to delete the selection from feature class "B" and the Append tool is used to append the selection from feature class "A" to feature class "B". In between, because we need a specific attribute of the features that is only contained in feature class "B", before we proceed with delete features from "B", we use Make Feature Layer to keep this selection and then Merge it with the selection coming from feature class "A" and Dissolve them on the unique field to keep one feature and the neccesary attributes.

The problem: When we run the model on a local geodatabase with 100 features (for testing) everything works fine. When we run the model on our production feature classes (Oracle SDE) with 500.000 features, the model crashes ArcMap when it reaches the Merge tool:

  "Error:  ArcMap.exe has encountered an error and needs to close. We apologize for the inconvenience" with the option to send an error report to ESRI.


texts: How to create dissolved buffer layer with pyqgis?

texts:
I have a line vector layer with a field name of "buffer_distance". I want create a dissolved buffer layer based on it's field. I write below code:

from qgis.utils import iface
from qgis.analysis import QgsGeometryAnalyzer
mc = iface.mapCanvas()
layer = mc.currentLayer()
Distance_Field = layer.fieldNameIndex("buffer_distance")
QgsGeometryAnalyzer().buffer(layer, "/tmp/buffer.shp", False, True, Distance_Field)

when run this code in QGIS 2.18 python console, It creates a new shapefile with all table fields of input layer but there is not any feature. What is problem?


texts: Create multiple pizza slice type shape from excel

texts: I am trying to make pizza type shape with various radious of each slice. Number of slices depend on number of years and radius of each slice depends on number of population. All I know the center longitude-latitude, and radius (num of population). Can anyone let me how to make such shape and how to design the excel sheet?


texts: Silence qgis.core error messages in reticulate

texts:
I am using the qgis.core Python module in R using reticulate. I start the QgsApplication using the following code.

# Start reticulate
reticulate::use_python(Sys.getenv("PYTHONHOME"), require = TRUE)

# Import modules
qgis.core <- reticulate::import("qgis.core")

# Create QGIS application
qgs <- qgis.core$QgsApplication(integer(), TRUE)

# Initiate QGIS application.
qgs$initQgis()

For whatever reason, I get the following messages when I run the qgs$initQgis() line:

ERROR: Opening of authentication db FAILED
ERROR: Unable to establish authentication database connection
ERROR: Auth db could not be created and opened

These are just printed to the R console, they do not interrupt the rest of the function or impact the rest of the code. These errors have also been reported as being inconsequential by other users: https://stackoverflow.com/a/38175594/5106834

It would be great if the error messages could not appear, but failing that, I would at least like to suppress the messages. So far I have tried suppressMessages(), sink(), invisible() and capture.output() and none have worked.

How can I silence these error messages?


texts: How to connect door nodes using a polyline

texts:
I am working on an indoor positioning project which need to transform CAD drawings to GIS shape files. I have transformed the file from CAD to ArcMap 10.1 already.

My current problem is that how could I transform the doors in CAD to a set of polyline to as an independent layer. Currently, a door is consisted of a arc and a line but the connected line (door sill) is actually what I want.

I want to know is there any simple geo-spatial method to do that. Hopefully it will not use programming because of limited deadline. I am sorry I cannot upload the picture at this stage.


texts: Copying/Projecting featureclass into a database and featuredataset of the same name fails

texts:
I've written a script in Python for ArcGIS 10.1 in order to migrate from an outdated coordinate system to a newer one while at the same time updating from personal geodatabases (.mdb) to file geodatabases (.gdb)- Everything works fine except:

When the tool tries to copy OR reproject (which I guess has to copy data as well) from a database containing a feature class with the same name as the featuredataset it fails, giving a 000171 Error 'Failed to open output workspace.

Changing the name of the feature dataset to something else will make it work (i.e. database.gdb\different\same
Changing the name of the feature class to something else will make it work (i.e. database.gdb\same\different)

I can see an easy workaround by just renaming stuff a couple of times, but I want to get to the bottom of it if possible- It's a very fringe case that doesn't affect many of the databases so I want to avoid more code just to handle this case if possible. It seems OK to create these duplicate names using normal methods (i.e. ArcCatalog won't stop you) and the inbuilt geo-processing tools in ArcToolbox (i.e. Project) don't fail in the same situation....it's just the arcpy library it seems?

(edit code added - Thanks to KimOs suggestion, the copy portion under the else statement now works- The projection still doesn't. All data has its projection defined correctly.)

            FCs =arcpy.ListFeatureClasses('*','All',dataset) # list feature classes within this feature dataset
            arcpy.AddMessage("Found these feature classes: " + str(FCs))
            for FC in FCs: # for the feature classes within a feature dataset
                descFC= arcpy.Describe(FC)
                if descFC.spatialReference.name == 'GD_1949_New_Zealand_Map_Grid': # reproject only NZMG data
                    arcpy.AddMessage(FC + " is NZMG, reprojecting to: " + newFD + descFC.baseName)
                    logging.info(FC + " is NZMG, reprojecting to: " + newFD + descFC.baseName)
                    arcpy.Project_management(FC, newFD + descFC.baseName, prjFile, 'New_Zealand_1949_To_NZGD_2000_3_NTv2', descFC.spatialReference)
                else: # else just copy it across
                    arcpy.AddMessage(FC + " is in " + str(descFC.spatialReference.name) + " , copying instead: " + newFD+ descFC.baseName)
                    logging.info(FC + " is in " + str(descFC.spatialReference.name) + " , copying instead: " + newFD+ descFC.baseName)
                    arcpy.Copy_management(FC,newFD + descFC.baseName, "FeatureClass")


texts: d.rast.edit error when attempting to edit a raster

texts:
I always get an error when I try to edit cell values in a raster with the GrassGIS command „d.rast.edit“.
The error under Windows XP/GrassGIS 6.4.3 says:

Exception in thread Thread-29:
Traceback (most recent call last):
  File "D:\OSGeo4W\apps\Python27\lib\threading.py", line
552, in __bootstrap_inner
    self.run()
  File "D:\OSGeo4W\apps\grass\grass-6.4.3RC1\etc\wxpython\co
re\gcmd.py", line 538, in run
    os.pathsep + env['PATH']
TypeError: 'NoneType' object is not subscriptable
Traceback (most recent call last):
  File "D:\OSGeo4W\apps\grass\grass-6.4.3RC1\etc\wxpython\gu
i_core\forms.py", line 657, in OnRun

gcmd.Command(cmd)
  File "D:\OSGeo4W\apps\grass\grass-6.4.3RC1\etc\wxpython\co
re\gcmd.py", line 389, in __init__

_("Error: ") + self.__GetError()))
core.gcmd
.
GException
Exception in thread Thread-30:
Traceback (most recent call last):
  File "D:\OSGeo4W\apps\Python27\lib\threading.py", line
552, in __bootstrap_inner
    self.run()
  File "D:\OSGeo4W\apps\grass\grass-6.4.3RC1\etc\wxpython\co
re\gcmd.py", line 538, in run
    os.pathsep + env['PATH']
TypeError: 'NoneType' object is not subscriptable
Traceback (most recent call last):
  File "D:\OSGeo4W\apps\grass\grass-6.4.3RC1\etc\wxpython\gu
i_core\forms.py", line 657, in OnRun

gcmd.Command(cmd)
  File "D:\OSGeo4W\apps\grass\grass-6.4.3RC1\etc\wxpython\co
re\gcmd.py", line 389, in __init__

_("Error: ") + self.__GetError()))
core.gcmd
.
GException

… and under openSUSE 11 SP2/ Grass GIS 6.4.2:

Traceback (most recent call last):
  File "/usr/local/grass-6.4.2/etc/wxpython/gui_modules/menu
form.py", line 662, in OnRun

gcmd.Command(cmd)
  File
"/usr/local/grass-6.4.2/etc/wxpython/gui_modules/gcmd.py",
line 350, in __init__

_("Error: ") + self.__GetError()))
gui_modules.gcmd
.
GException
:
Execution failed: 'd.rast.edit --overwrite --verbose
input=DEM_m_float_polypoi@PERMANENT output=123456'
Details:
Error:
Traceback (most recent call last):
  File "/usr/local/grass-6.4.2/etc/wxpython/gui_modules/menu
form.py", line 662, in OnRun

gcmd.Command(cmd)
  File
"/usr/local/grass-6.4.2/etc/wxpython/gui_modules/gcmd.py",
line 350, in __init__

_("Error: ") + self.__GetError()))
gui_modules.gcmd
.
GException
:
Execution failed: 'd.rast.edit --overwrite --verbose
input=DEM_m_float_polypoi@PERMANENT output=123456'
Details:
Error:

Does anyone have an idea about the reason why I'm getting this error and/or how to avoid it ?


texts: PostGIS query losing geometry column information

texts:
Using PostGIS I want to collect all attributes into one polygon. The old table looks like this:

Column   | Type
---------+------
xyz_cols | `xyz`
geometry | geometry(MultipolgyonZ, 2193)

Performing the query:

CREATE TABLE new_table AS
SELECT
  1::INT as id,
  ST_Collect(geom)
FROM
  old_table

which works but I lose all geometry information and the output table is like this:

Column   | Type
---------+------
id       | INT
geometry | geometry

Even if I wrap the ST_Collect in a ST_Multi() I lose the geometry information.


texts:
I am using ArcGIS 10 and have a polygon shapefile in which many of the numeric fields are Null.

The polygon shapefile in question is produced by joining a polygon shapefile with a point shapefile and the "Count_" field contains the count of how many points fall inside the polygon. I would like to convert these Nulls to 0. So I have used two scripts in Field Calculator(Python). Both of the scripts execute successfully but produce no results i.e. the null values are not replaced. The two scripts are as follows:

def rep_null(count):
  if count is None:
    return 0
  else:
    return count

and

def rep_null(count):
  out=0
  if count>=1:
    out=count
  return out

What is the way out?


texts: ArcMap - Driving distance between 2 sets of points

texts:
I have a question regarding the calculating driving distances between any two points.

I have 2 CSV files, an origin and a destination file, which look like the following:

Origin                                    Destination
ID       Latitude     Longitude           ID       Latitude     Longitude
1        33.9425      -118.408,1          15       44.08        -69.05
2        66.0581      -23.1353,2          16       46.41        -68.49
                                          17       44.26        -70.49

The goal is to calculate the driving distance between any origin and each destination, i.e., given the two CSV files, the result should contain 6 routes. While I could solve this problem using Google Maps or Mapquest, I tried to go with ArcMap, because the original origin CSV file contains 500 lines, whil the destination file contains around 42.000 lines, building around 20 million routes.

When opening a new project, I choose the "Coterminous United States" map as a map, because all points in both files are located in the United States.

Could it be possible that there are no underlying routes for the Network Analyst to take into account?

If so, how can I generate them?

Importing the data sets worked fine (via File -> Add Data -> Add Data). However, displaying them on the map showed, that there was no object id - so I exported them and re - displayed them on the map, leading to the following points:

According to this GeoNet answer, adding a route should be possible now. However, after activating the Network Analyst, only the "Network Analyst Options" are clickable.

How can I tell the Network Analyst to calculate the distances between both sets?


texts: Overpass API call JSON result is different from XML result

texts:
I'm calling this:
http://www.overpass-api.de/api/xapi?debug=*[amenity=hospital][bbox=13.20524,43.70861,13.22842,43.72338]

And I have a result that is correct, since it shows me the only hospital in the area. Then, following the instructions explained in this answer, I've modified the call like this:
https://www.overpass-api.de/api/interpreter?data=[out:json];node[amenity=hospital](43.70861,13.20524,43.72338,13.22842);out%20meta;

With this I do not have the same result, since "elements" is empty, like there is no hospital in the area. If I remove the "amenity=hospital" filter, I do have correct nodes from the area selected (but still, not the hospital!), so why the results are different?
Also, If I wanted to use the first solution, how could I read the xapi file I get? (I'm using angular for a web app)


texts: sampling error using gdalwarp on a subsetted VRT?

texts:
I am currently trying to prepare an SRTM subset including height conversion from EGM96 to WGS84 using gdalwarp. I encountered some strange shifts in part of the image.
Here is a reproducible example:
Prepare the necessary file for the height conversion:

cd /usr/share/proj
sudo wget https://download.osgeo.org/proj/vdatum/egm96_15/egm96_15.gtx
sudo chmod 644 egm96_15.gtx

Download two SRTM 3Sec tiles and combine them in a VRT:

wget http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_39_03.zip
wget http://srtm.csi.cgiar.org/wp-content/uploads/files/srtm_5x5/TIFF/srtm_38_03.zip
gdalbuildvrt dem.vrt /vsizip/srtm_39_03.zip/srtm_39_03.tif \
/vsizip/srtm_38_03.zip/srtm_38_03.tif \
-te 9.497235 45.432812 13.183811 47.336189

Use gdalwarp to create two new files, one that is just a simple conversion of the original file content and one that contains heights converted to WGS84 ellipsoid:

gdalwarp -r bilinear -of GTiff dem.vrt dem.tif
gdalwarp -r bilinear -of GTiff -s_srs EPSG:4326+5773 -t_srs EPSG:4326 \
dem.vrt dem_ellp.tif

Compute their difference as GTiff and JPEG to get the EGM96-WGS84 difference:

gdal_calc.py -A dem.tif -B dem_ellp.tif --outfile=dem_ellp_diff.tif --calc="A-B"
gdal_translate -of JPEG -scale -outsize 50% 50% \
dem_ellp_diff.tif dem_ellp_diff.jpg

The output looks like this:

To the West and East the image looks like expected with simple gradient lines. In the center however it looks like a DEM aspect file. On closer inspection of the two files, one can see a pixel shift to the East that is introduced by gdalwarp in the second (ellipsoid height) file.
Is this problem known to someone? It does not appear when the DEM is not subsetted during creating the vrt file, i.e. calling gdalbuildvrt without the te argument:

Also, the images look the same when using a different spatial CRS e.g. UTM
via
-t_srs EPSG:32632

The problem is also not present when alternatively setting the extent via gdalwarp instead of gdalbuildvrt:

gdalbuildvrt dem2.vrt /vsizip/srtm_39_03.zip/srtm_39_03.tif \
/vsizip/srtm_38_03.zip/srtm_38_03.tif
gdalwarp -r bilinear -of GTiff -te 9.497235 45.432812 13.183811 47.336189 \
dem2.vrt dem2.tif
gdalwarp -r bilinear -of GTiff -s_srs EPSG:4326+5773 -t_srs EPSG:4326 \
-te 9.497235 45.432812 13.183811 47.336189 dem2.vrt dem2_ellp.tif
gdal_calc.py -A dem2.tif -B dem2_ellp.tif --outfile=dem2_ellp_diff.tif \
--calc="A-B"
gdal_translate -of JPEG -scale -outsize 50% 50% \
dem2_ellp_diff.tif dem2_ellp_diff.jpg

GDAL v2.2.3, Ubuntu 18.04, libproj v4.9.3-2


texts: Getting a 3D view from Google Maps into QGIS?

texts: I guess it would be safe to assume that there is no straight-forward way of doing this, since 3D functionality was just rolled into QGIS, but is there a way of importing the 3D (tilted) view from the online google maps as a dynamic 3D basemap into QGIS?


texts:
I've found a few questions similar but think mine is much more simple.
I would like to run: Export -> Save As as a batch command.
I have 1 layer that has 9595 DEM raster overlaid and I would like to save the 1 layer underneath, clipped into the extents of each DEM file.
Is this possible?


texts: How To Loop in A Dataset and Delete Specific Fields in All Feature Classes With ArcPY

texts: How can I use ArcPy to loop into all Feature Classes in a dataset Structurein  Geodatabase GEWA.gdb and delete Two fields of them called MapBy and EnabledBy if they exist? (All of feature classes have the Mapby but only some of them has the EnabledBy).


texts: Esri Leaflet identify() only on added layer from MapService

texts:
I try to identify features only from the wastepit layer. In this example I am adding two layers: wells and wastepit. I am only calling the identify function on the wastepit layer but I still get features if I press on a location for a well? Please see the codepen for live example. How can I do so I am only get the featureInfo from the wastepit layer? I dont understand why I in fact get feature from all (8) layers in the MapService?

    wells = L.esri.dynamicMapLayer({
            url: url,
            layers: [0]
          }).addTo(map);

pipe = L.esri.dynamicMapLayer({
                url: url,
                layers: [4]
              }).addTo(map);

wastepit = L.esri.dynamicMapLayer({
            url: url,
            layers: [5]
          }).addTo(map);

function getFeatureInfo(e) {
  wastepit.identify()
    .on(map)
    .at(e.latlng)
    .run(function(error, featureCollection){

https://codepen.io/GitMorin/pen/MLQvdM?editors=0010


texts: Are there any web services providing country border data?

texts:
I am an Android developer (working in Java) and want to write a geographical quiz game. The quiz will need a list of countries and will randomly choose one to highlight on a map (probably using Google Maps). Then the user can type in the name of the country.

What web services are available that will give me the polygon data of the border of a given country?


texts: How to automatically run some scripts for QGIS python console by defining the python macros checkbox

texts:
I'm working with QGIS 2.8 in windows 10.
I've defined the python macros checkbox to QGIS start and simultaneously run a script that has the same path from my "QGIS.qgs" file.
In the Project>Project properties>Macros I defined:

import myScript1

def openProject():
    myScript1.run()

def saveProject():
    pass

def closeProject():
    pass

Then, I've wrote a script (myScript1.py) in the same path from the QGIS project.

That worked as I wanted. Now I want to do the same for another script that should run in a sequence of the first one. Imagine I wanted to run first the "myScript1.py", then "myScript2.py". How should I define the python macros checkbox from the QGIS project? Should I write the "myScript2.py" in the same path of "myScript1.py"? The "myScript1.py" is to load a raster and define the legend:

from qgis.core import *
from PyQt4.QtGui import *
from PyQt4 import QtGui
from PyQt4.QtCore import *
from PyQt4.QtXml import *
from PyQt4.QtCore import QSize
from PyQt4 import QtCore
from qgis import core, gui
from qgis.utils import iface

def run():
    lyr = QgsRasterLayer("C:/PATRICIA/IG/TESTE/TEST/FWI_20150802.dat", "FWI_20150802")
    s = QgsRasterShader()
    c = QgsColorRampShader()
    c.setColorRampType(QgsColorRampShader.INTERPOLATED)
    i = []

    i.append(QgsColorRampShader.ColorRampItem(8.5, QtGui.QColor('#2b9e27'), '8.5'))
    i.append(QgsColorRampShader.ColorRampItem(17.2, QtGui.QColor('#fee717'), '17.2'))
    i.append(QgsColorRampShader.ColorRampItem(24.6, QtGui.QColor('#ec7000'), '24.6'))
    i.append(QgsColorRampShader.ColorRampItem(38.3, QtGui.QColor('#cc292c'), '38.3'))
    i.append(QgsColorRampShader.ColorRampItem(100, QtGui.QColor('#75263b'), '100'))

    c.setColorRampItemList(i)
    s.setRasterShaderFunction(c)
    ps = QgsSingleBandPseudoColorRenderer(lyr.dataProvider(), 1, s)
    lyr.setRenderer(ps)
    QgsMapLayerRegistry.instance().addMapLayer(lyr)

The next script that I want to run in sequence of the above one is "myScript2.py" and is to create a print composer based in a defined "template.qpt" and convert it to an image output (png):

from qgis.core import *
from PyQt4.QtGui import *
from PyQt4 import QtGui
from PyQt4.QtCore import *
from PyQt4.QtXml import *
from PyQt4.QtCore import QSize
from PyQt4 import QtCore
from qgis import core, gui
from qgis.utils import iface

def run():
    layers =iface.legendInterface().layers()
    canvas=iface.mapCanvas()

    for layer in layers:
        iface.legendInterface().setLayerVisible(layer, True)
        myFile = "C:/PATRICIA/IG/TESTE/Template.qpt"
        myTemplateFile = file(myFile, 'rt')
        myTemplateContent = myTemplateFile.read()
        myTemplateFile.close()
        myDocument = QDomDocument()
        myDocument.setContent(myTemplateContent, False)
        newcomp = iface.createNewComposer()
        newcomp.composition().loadFromTemplate(myDocument)
        newcomp.composition().refreshItems()

        c = iface.activeComposers()[0].composition()
        image = c.printPageAsRaster(0)
        image.save('C:/PATRICIA/IG/TESTE/OUTPUT/FWI_Arome.png','png')

Does anyone can help me?


texts: How to (image) capture map and copy it to clipboard using QGIS?

texts: Essentially, I want to copy a map to clipboard (not Save as Image).


texts: Raster algebra to tiled maps on a sever with open source solutions

texts:
Let's say I have a collection of data in raster format.  It changes daily.  Once each day I would like to do some fairly complex raster algebra on these datasets.  I would like to be able to do the calculations automatically (so I don't have to manually run the calculations).  Ideal I can use python to code the calculations.  I would like to take the results of the raster calculations and then tile them up to be used as a service displayed on my website.  Again, I'd like to do this automatically, so that I don't have to set up the tiles each day.  I am pretty sure I can do this with some collection ESRI products, but it would be great to have an open solution.  I suppose that once I have the results of the raster algebra, I could use geoserver or similar to distribute the results, but a more streamlined solution would also be good.

Bottom Line: What software can I use to do raster algebra calculations on a sever?


texts: Performing a spatial join to obtain attribute field summaries of nested polygons that are located within polygon regions using ArcPy

texts:
I have created an ArcPy script tool that finds depressions from an input DEM, and then generates closed contours (both polyline and polygon outputs) as well as regions (raster and polygon outputs).  I need my tool to be able to output a final layer that provides summary attributes of the closed contour polygons that fall within each region from the regions polygon feature class.
To do this, I know I need to perform a spatial join with the closed contour (cc) polygons and the regions polygons to obtain statistics from the attributes of the cc polygons within the regions; however, I want multiple statistics for each attribute field, such as 'Max','Min',and 'Range'.  I can successfully do this using the 'Spatial Join' within ArcMap 10.6's user interface, but I want to do this using ArcPy - I think I need to use fieldmappings to do this and set up merge rules, but I'm not sure how to do this step.  I have used fieldmappings to add both tables, and to remove fields that I don't want to retain, but I am not sure where to proceed from there.
I have included a code snippet of my current fieldmappings before performing the spatial join below:
outFeatures = os.path.join(outLoc,outElim+ "_Final")
fieldmappings = arcpy.FieldMappings()

fieldmappings.addTable(outRegions)
fieldmappings.addTable(outElim)

keepers = ["area_sqft","perim_ft","Min_Z","Max_Z","Range_Z","circularityIndex","solidity","gridcode","length_ft","MBG_Width","MBG_Length","MBG_Orientation"]

#remove all output fields that aren't desired, then perform spatial joins
for field in fieldmappings.fields:
        if field.name not in keepers:
                fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex(field.name))

#set merge rules here to obtain min, max, and range for kept attribute fields

arcpy.SpatialJoin_analysis(outRegions,outElim,outFeatures,"JOIN_ONE_TO_ONE","KEEP_COMMON",fieldmappings,"CONTAINS","","")


texts: Getting Graphic from Graphics Layer using ArcGIS API for JavaScript

texts:
How do I get a Graphic from the Graphic Layer like the way one can implement in OpenLayers using getFeatureById()?
Presently, I am attaching attributes to graphic and then by getting graphic layer by layerId, I loop through all the graphics in that layer until the loop satisfy the condition and returns me the correct Graphic I am looking for.
So it works in following way,
`//This returns graphics layer from the list of graphics layer
var graphicsLayer = map.getLayer(layerId);
//Not supported by API
//var graphic = graphicsLayer.getGraphic(graphicId);`

`//Current approach
for(i=0;i<graphisLayer.graphics.length;i++){
    if(graphisLayer.graphics[i].attribute.fid == myFid){
        //Do something with this graphic
        break;
    }
}`

However, I am seeking another approach.


texts: Invalid filter when using MapServer as a WFS

texts:
I'm trying to use Openlayers and Mapserver to display a shapefile. My HTML file is as follows:

  <title>OpenLayers Example</title>

This loads a google-map base map, and should overlay a shapefile over all of England. When the page loads the coordinates of the bounding box of the map are sent to MapServer. I can see that it is sending just the lon/lat values. In my mapfile, I have the projection of the data defined as follows, in both the intro part, and also in the layer:

 PROJECTION
    "init=epsg:27700"
 END

However, the response from the server gives me an error WFS server error. Invalid or Unsupported FILTER in
GetFeature, and I don't know whether it is due to an incorrect bounding box, or whether I have a mistake with my query using PropertyName. From searching online this appears to have been a bug, but I don't know whether it is still a problem or not.

<?xml version="1.0" encoding="ISO-8859-1"?>
<ows:ExceptionReport xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xmlns:ows="http://www.opengis.net/ows" version="1.1.0" language="en-US"
xsi:schemaLocation="http://www.opengis.net/ows http://schemas.opengis.net/ows/1.0.0/owsExceptionReport.xsd">
  <ows:Exception exceptionCode="InvalidParameterValue" locator="filter">
    <ows:ExceptionText>msWFSGetFeature(): WFS server error. Invalid or Unsupported FILTER in
GetFeature : &lt;
    Filter ogc="http://www.opengis.net/ogc"&gt;
  &lt;BBOX&gt;
    &lt;Envelope gml="http://www.opengis.net/gml" srsName="EPSG:4326"&gt;
      &lt;lowerCorner&gt;-55.50109375 43.134977942511&lt;/lowerCorner&gt;
      &lt;upperCorner&gt;55.24109375 59.593781185019&lt;/upperCorner&gt;
    &lt;/Envelope&gt;
  &lt;/BBOX&gt;
&lt;/Filter&gt;
</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>

What do I need to do to fix this? Should I make sure that I am passing a bounding box in epsg:27700 coordinates? Or can Mapserver accept lon/lat and project them into the required coordinate system?

I have asked another question on the same topic recently, but I had a different problem then!


texts: Input parameters not saved in HEC-RAS 5.0.6

texts:
I have worked in HEC-RAS before, even in the version I am trying to use right now (5.0.6) but I have never encountered this issue. It seems I am unable to introduce decimal numbers in HEC-RAS, due to the fact that the program saves float numbers as integers. I have checked the Regional Settings in the system, I've changed the measure unit (just to check if that is the issue) but without success.

Also, this symbology error come up along with the aforementioned issue:

Has someone experienced this issue?

Later edit: I think I've found the solution. The issue what that, in the Regional Settings, the decimal separator was a comma and not a dot. Now it seems to work well.


texts:
I need to do some raster calculation, I've been trying to make it work but I am not sure everything is working properly.

I need to apply a formula, but I am unsure if I am doing it properly. The formula is needed to find out the luminance:

  L(Bx) = (CN(Bx) – 1) * G(Bx)   -  Q0(BX)

Bx is the band, G is 0.64 and Qo is -5, CN is the number of the band ( in this case 3), the name of the file is "B30@1"

Should it look like this ?

  ( 3B30@1 - 1 )  * 0.64B30@1 - -5B30@1

After I apply this formula, the image become a gray square, is this normal? Have I applied the formula correctly


texts: configuring cartodb for running on local machine

texts: How to  configure cartodb to run fully on local machine. I have setup a development version of cartodb. In this the tables are limited to 5, the disk space is limited and the view. Is these variables can be configured to be run on local server. I have looked the database.yml and app_config.yml. Is there any docs available for cartodb server configuration on  local machine.


texts: How to set interactivity on multiple columns in CartoDB?

texts:
I'm trying to create a sublayer with interactivity on it for a tooltip. The tooltip is going to need data from a column named plant_name. So I am setting the interactivity property as I create the sublayer, like this:

        sublayer = layer.createSubLayer(
          sql: sql,
          cartocss: css
          interactivity: ["cartodb_id", "plant_name"]
        )

The problem is I get an error saying "layer1: Postgis Plugin: ERROR: column "plant_name" does not exist." even though the column is in the table.

What am I doing wrong?


texts: Minimum and maximum values when raster statistics are lacking

texts:
How does ArcMap determine minimum and maximum raster values (e.g. for rendering colour ramps) when a raster file lacks an auxiliary file (e.g. aux.xml) describing statistics?

This question is related to my previous question here, which describes a situation where, for GeoTiffs generated by gdalwarp, ArcMap displayed the minimum as 268.048, while the true minimum was 91.2931.

Does ArcMap inspect an initial subset of cells, and set statistics based on those cells' values? If so, can this behaviour be prevented other than by ensuring auxiliary files are created and distributed with the raster? It seems like risky behaviour, since unless a user ensures that statistics have been fully calculated, rendered output might easily be misinterpreted.


texts: Responsive d3 choropleth showing commuter flows between counties?

texts:
I would like to create a responsive choropleth using d3. When one clicks on a county, the choropleth updates to reflect commuter flows between that county and all others. Example input data looks like this:

"home","work","flow"
01001,01007,10
01001,01015,10
01001,01021,383
01001,01047,271
01001,01051,1624
01001,01073,71
01001,01081,50
01001,01097,80
01001,37119,10
01001,40143,10
01001,48029,45
01001,48141,10
01001,51041,35
13121,01003,35
13121,01071,15
13121,01073,69
13121,01081,4
13121,01089,70
13121,01097,55
13121,01121,10
13121,01125,15
13121,04005,10
13121,04013,65

In the unemployment choropleth example where the choropleth map is static, the input data is read in like this:

queue()
    .defer(d3.json, "data/us.json")
    .defer(d3.tsv, "data/unemployment.tsv", function(d) { rateById.set(d.id, +d.rate); })
    .await(ready);

Future calls on rateById use get() (e.g., rateById.get(d.id)). See below from the unemployment example.

quantize(rateById.get(d.id));

If I have the commuter flow two-dimensional dataset, I expect that calls on flowByHomeByWork would look like flowByHomeByWork.get(currFIPS).get(d.id). For example,

quantize(flowByHomeByWork.get(currFIPS).get(d.id));

How do I map the data when reading it in with queue() to facilitate this?


texts: Wedge neighborhood on Focal Statistics calculation

texts: I am using Focal Statistics in Spatial Analyst of ArcMap10.2, and use the wedge type of neighborhood for a radius of 10 km around each processing pixel. I want to calculate the maximum elevation for four sectors of aspect which are: Northeast, Southeast, Southwest, and Northwest, which in degrees would mean 0-90, 90-180, 180-270, 270-360 respectively, If 0 degrees were on the north and the wedge would extend clockwise. Of course, as I have read from the ArcGIS Help, "The wedge extends counterclockwise from the starting angle to the ending angle. Angles are specified in arithmetic degrees from 0 to 360, where 0 is on the positive x-axis (3:00 on a clock)". Nevertheless, what I have observed is that 0 degrees are placed on the west (9:00 on a clock) and the wedge extends counterclockwise (90 degrees placed on 6:00 o clock, 180 placed on 3:00 o clock, 270 placed on 12:00 o clock). To be more precise, the images below depict a part of my region. The first image has a wedge with start angle at 0 degrees and end angle at 90, the second has an angle of 90-180, the third 180-270, and the fourth image has an angle of 270-360 degrees.


texts: Setting Data Driven Pages PDF File Name?

texts:
I have a mapbook I created using ArcGIS Desktop 10.0 that has its index layer as a 4 digit number.

I just wanted to export my maps as PDFs so they only show the 4 digit number, at the moment they are showing the name of my mxd and the 4 digit number.

e.g. MyMapbook_5132

Is there a simply way I can export them to only show the numbers,

e.g. 5132.pdf rather than MyMapbook_5132.pdf

PS - my 4 digit number comes from a field within my index grid


texts: What can I do with a read-only Feature Service that I cannot do with a Map Service?

texts:
If I have a feature service with the editing capabilities disabled what can I do with it that I cannot accomplish with a regular Map Service?

I can't think of a purpose for having a Feature Service except for editing.  Even if I was building a JS API application and I wanted the geometries I could create Feature Layer object using a MapService URL which would query the service to get the geometries.  On top of that with the MapService I have the option of just using a regular dynamic layer in which case I get only the image back.

I also tried using the read-only feature service as input to a geoprocessing tool and that doesn't work either.

The regular Map Service seems to be a lot more flexible than the Feature Service for everything other than editing.


texts: How to check if my OpenLayers Map has finished to render?

texts:
I am exporting an image from the current view of my map, but sometimes it exports an incomplete image like this:

And I need the image looks like this:

This second image I took a screen shot of my map to illustrate what I need.
I'm using OpenLayers.


texts:
I was trying to filter openstreet map data using osmosis. I used the following

osmosis --read-xml city.osm --tf accept-ways highway=* --used-node --write-xml highways.osm

But while i run this command in osmosis with osm file "city.osm" it throws run time exception

SEVERE: Thread for task 1-read-xml failed org.openstreetmap.osmosis.core.OsmosisRuntimeException: Cannot represent 2169105896 as an integer. at org.openstreetmap.osmosis.core.util.LongAsInt.longToInt(LongAsInt.java:33)

What could be the problem?

I updated osmosis to 0.43.1 version and it shows runtime exception

    SEVERE: Thread for task 1-read-xml failed
    org.openstreetmap.osmosis.core.OsmosisRuntimeException: The entity timestamp attribute is missing.
at org.openstreetmap.osmosis.core.domain.common.UnparsedTimestampContainer<init>(UnparsedTimestampContainer.java:34)


texts: Merging Feature Classes and keeping the metadata

texts:
Is there a way that while merging two feature classes we can keep the metadata?

We are merging feature classes with templates so we can import specific fields into our data.

Is there a way to keep a metadata template that is imported into our feature class templates with our existing feature class?


texts: Blue indentation error in QGIS' Python Console

texts:
I have small problem with Python Console in QGIS.
Normally, I write some code in text editior and then copy whole code to python console. Unfortunately, the console underlines the indentations with blue lines and I got en error: "unindent does not match any outer indentation level."
To solve it I need to copy paste tabs from lines where it is ok.
My question is: why is this problem occuring and how to solve in easier way? I need to say that copying tabs is unnecessary waste of time.
I am attaching the picture to illustrate the problem.


texts: Improving poor performance of ArcPy update cursor with ArcSDE and Editor?

texts:
I am using ArcGIS 10.5 in an ArcSDE environment.

I have two tables, they are the exact same table schemas, but one is being edited and those edits are being pushed to the main table.  I am a little weak in my python (but learning) and I just started to work with da.cursors and I have gotten this script to work,  but both tables have 3800 records in them,  and when this runs to update it takes 4 hours.  This is just way too long imho.

I am pretty sure I am doing something wrong here but I dont know what it is.  Also I am open to any other suggestions on how to best accomplish this task.  Some people have suggested joining the two tables based on the UUID and then using the calculate tool to copy over all values,  some have suggested using a dictionary.

I also think I dont need the whole arcpy.da.Editor part, but that was in the code sample I found, so I was a little confused.

Please dont feel like I am asking anyone to do this for me,  I have gotten this to actually work, but it works very poorly and I am afraid my logic is not very elegant and may actually cause problems later on (like locked dbs, etc...)  I would just like a critique of the code and any suggestions on how to better accomplish this.

here is my code

Dump_FC = r'C:\Users\ME\Documents\ArcGIS\Default.gdb\CondemnationTimeline_UUID_xls_CondemnationTimeline_UUID'
Main_FC =  r'Database Connections\IS memvmtestsql3.XXX.ad.sde\IS_GDB.GISADMIN.PW_CodeEnforcement_CondemnationTimeline'

#Load all field names into list named field_names (crazy huh?)
dsc = arcpy.Describe(Dump_FC)
fields = dsc.fields

field_names = []

#Remove any field named ObjectID or ObjectID_1 and create list of fields

for field in fields:
    if field.name != "OBJECTID" and field.name != "OBJECTID_1" :
        print field.name + " != " + dsc.OIDFieldName + " " + str(field.name != dsc.OIDFieldName)
        field_names.append(field.name)
    else:
        continue

fc1 = r'C:\Users\XXX\Documents\ArcGIS\Default.gdb\CondemnationTimeline_UUID_xls_CondemnationTimeline_UUID'
fc1_fields = field_names
fc2 = r'Database Connections\IS memvmtestsql3.XXXX.ad.sde\IS_GDB.GISADMIN.PW_CodeEnforcement_CondemnationTimeline'
fc2_fields = field_names

# The Editor class allows use of edit sessions and operations to manage database transactions.
edit = arcpy.da.Editor(workspace)
edit.startEditing(True, False)
edit.startOperation()

i=1

print "starting main table update"

# Open an da.updateCursor on FC1 starting at ROW1 , at ROW1 open second da.updateCursor on FC2 and go ROW by ROW looking for matching GRID #

with arcpy.da.UpdateCursor(fc1, fc1_fields) as cur:
    for row in cur:
        cur2 = arcpy.da.UpdateCursor(fc2, fc2_fields)
        print "on row ==================== " + str(i) + "   " + str(row[64])
        i = i + 1
        for row2 in cur2:
            print row2[64]
            if row[64] == row2[64]:   # if UUIDS in FC1 match FC2 then copy over all columns
                print "*** match found ***"

                j = 0
                while j < len(field_names):
                    row2[j] = row[j]
                    j = j + 1

            cur2.updateRow(row2)

del cur, cur2, row, row2

edit.stopOperation()
edit.stopEditing(True)


texts: Checking if attribute field for labeling is present to avoid ArcMap Drawing Errors?

texts:
I have a project with several layers in the TOC that point to feature classes or shapefiles. Each of them is using labels and everything is ok.
If I replace the source of these layers with some other feature classes/shapefiles with different attribute fields, sometimes it happens that the atribute field for labeling is in the new source data missing. That means that in the layer properties the Label features in this layer checkbox is ticked but the window with the attribute field is blank. Then when refreshing the layout view, the ArcMap sometimes (I don´t know why only sometimes but still) throws an error window like this:

Since I can´t rely on ArcMap because the error window sometimes pops up and sometimes not, is there a way to use arcpy for looping through all layers in the TOC and find those with broken labels (=missing attribute field for labeling)? I can´t find anything about this possibility and no one has asked similar question.


texts:
I am a QGIS newbie trying to work through Lex Berman's QGIS tutorial.  I'm stuck in the "join" unit, because I do not see a Vector menu on my menu bar.  The screen shots show that there should be a Vector menu after the Help menu and before the Window menu on the menu bar.

Is this a plugin I need to install, or did I install the software improperly?  (I've found that I've had to install a couple of plugins myself, things that aren't described outright in the tutorial.)


texts: Map disappearing when adding graphic layers using Web ADF of ArcGIS Server?

texts:
I'm trying to add a graphic layer in a web ADF application with arcGIS for asp.net, but every time I add a graphic layer to the MapResourceManager, the map disappears!

To better understand this: here's the simple app with the only element in the MapResourceManager being the map service (OGC (WMS) Service) which provides me the map

http://imageshack.us/f/715/pro2w.jpg/

And here's the MapResourceManager with another graphic layer, the map misteriously disappears

http://imageshack.us/f/543/pro1y.jpg/

What's wrong with that?


texts:
Does anyone know how to fix this ?

I used to have QGIS 2.2 and 2.4. I've completely removed 2.4 and installed 2.6 by OSgeo4W (Advanced installation). Everything, like GRASS, tora, udig is working fine except QGIS. Moreover if i open 2.2, the plugins are not working. Might be due to python compatibility i guess.

I've tried the steps mentioned here How to fix QGIS error "Entry Point could not be located"?. But didn't work out.

Edit:

I've completely uninstalled OSgeo4W and again followed the same - Advanced install - option. This time the error was ippj-5.3.dll missing. Again i've uninstalled every OSgeo4W component and also thoroughly removed registry traces.

Now, i've installed using Express install option and it was installed successfully & QGIS is running properly.


texts: Getting info from multiple rows in single popup of CartoDB?

texts:
If you have a table like this:

Name    | Item   | State
---------|--------|-----------
 Alice   | Apple  | Alabama
 Ben     | Banana | Alabama
 Caitlin | Cherry | California

Assuming you have created polygons based on the states, how could you create an 'info window' which combines the details of all rows in the same state?

So for example, the popup for Alabama might be something like:

State: Alabama
People: 2
Items: Apple, Banana

As far as I can tell, the default behaviour for Cartodb would be to place one polygon on top of the other, and for the info window to display the info from the top most polygon only.


texts: select parcels in polygons and add fiels depending on the polygons

texts:
I have a lot of feature classes (results of service area with different polygons (0 to 10/ 10 to 20/ 20 to 40 and 40 to 60) (figure 1)

Also i have a shapefile with a lot of parcels (figure 2)

Now i want to select the parcels that are located in the polygons (figure 1)
thereafter i want to add a field in the shapefile 'parcels' with a number 4 if the parcel lies in the polygon of 0 to 10, a number 3 if the parcel lies in the polygon 10 to 20, number 2 for polygon 20 to 40 and number 1 for polygon 40 to 60.

For my next result of service area (figure 1), so from another timestep or location, i also want to select the parcels that are located in this polygons, add a field with number 4 if.., number 3 if...
So i received a table like figure 3


texts: Setting COLORRANGE for rasters with variable minimum and maximum value in MapServer?

texts:
I have a series of rasters picturing the same phenomenon/product eg. the accumulated precipitation and I use the same mapfile for all. Naturally the minimum is always zero but maximum is variable from dataset to dataset. But so far I know only how to set the COLORRANGE directive to fixed ranges of values like for example below:

CLASS
  EXPRESSION ([pixel] >= 15 AND [pixel] < 20)
  STYLE
    COLORRANGE 74 120 182 50 87 160
    DATARANGE 15 20
    RANGEITEM "pixel"
  END
END

So, is it possible to replace the 15 and 20 from the above example with the raster dependant minimum and maximum pixels value?

I know that I could probably set it by obtaining the values prior to the setting of the WMS source and sending them together with the other source parameters to the CGI as min=0... and setting them n the mapfile via %min%, %max%. But I am interested if mapfile is able to extract this info on its own.


texts: Free GeoSserver SLD for 2009 TIGER line data (preferably Google Maps color scheme)

texts:
I need to render a base map using OpenLayers. I have downloaded 2009 TIGER line data. I need to get up and running quickly. I think I can save some time if I use ready made SLDs, instead of designing them my elf.

I am looking for some decent looking SLDs especially for the road network and its labels. In fact I prefer if the SLDs are Google Maps style color scheme!

Please point me to links where I can download these SLDs.


texts: I have a set of points (about 400) shp file  which represents cities and the shp file of the road network. Also I have a new shp file of obstacles on the road. Now I want to count obstacles on the shortest routes between all of 400 places. Does somebody know how to do this?


texts: Using ModelBuilder Iterator sub-model then append?

texts:
I've created an iterator as a sub-model that collects values, then a main model that appends those values together upon completion of all iterations. My issue is if the iterator runs 10 times, then I get 10 appended records of the 10th iteration (the final value) not the 1 - 10 unique values. If I remove the Collect Values function and just add the append function there the model works perfectly, but that isn't suitable since I don't want to have to append each and every iteration to save time:

1) Iterator Sub Model

2) Main Model


texts: QGIS web client - Setting up WSGI search (mod_wsgi)

texts:
I am trying to set up the QGIS Web Clients WSGI search. I've succefully set up my QGIS WebClient on a windows OS. Now I tried installing the mod_wsgi.so module in the Apache modules directory, added the line LoadModule modules wsgi_module / mod_wsgi.so always in apache but as soon as I plug in the conf apache:

 #mod_wsgi

 WSGIScriptAlias / WSGI "C: \ OSGeo4W \ apache \ htdocs \ qgis-web-client-master \ WSGI"      WSGIScriptAliasMatch ^ / WSGI / ([^ /] +) "C: \ OSGeo4W \ apache \ htdocs \ qgis -Web-client-master \ $ 1.wsgi"

The apache server does not reboot. I'd like to understand what am I doing wrong and how to further set up my apache to execute wsgi files?


texts:
I have done that many times and it usually works but I'm facing a problem now that I do not understand. I'm trying to load a raster (.tif) with a EPSG4326 and it is a one band file in my database from the windows command prompt. I use the following function:

raster2pgsql -s 4326 -I -C -M -R climat.tif -F -t 25x25 clim.climat | psql -h localhost -U postgres -d biotic_interactions

it seems at first glance that it is working as I get 0 and 1 inserted and that my new table appears in pgAdmin. However, when I open it in QGIS and I confront the table that I just imported to the initial raster I see the value of the pixels are not the same than in the original file. Also, I noticed that while importing the raster, I get a few NOTICE:

  return error message: numeric field overflow
  unable to add constraint: nodata_values
  unable to add constraint: enforce_nodata_values_Rast

I'm using PostgreSQL 9.4.


texts: Calculating WGS84 (EPSG:4326) area in C# or SQLGeometry?

texts:
I am not normally a GIS developer, I am just helping out with a bug and I am stumped on what is going on here, could I get some help please.

I have two calculations, one in 4326 (WGS84) and one in 3857 (Google Mercator).  The coordinates were converted from 4326 -> 3857 and I have confirmed those are correct with ArcGIS and various online site conversion tools.  The results I have are...

-114.50,38.0=-12746081.696,4579425.813

-114.00,40.0=-12690421.950,4865942.280

If I do an estimate of the area for 3857 by calculating the difference in the x and y points then multiplying the differences and dividing by 1 mil to get sq km, I get the result just under 16,000 sq km.  Running SqlGeometry.STArea on a polygon bbox using these points this gives me a similar number.

If I convert the lat long to UTM I receive the numbers ...

-114.50,38.0=719510.335753775,4208764.46330096

-114.00,40.0=243900.352021924,4432069.05663229

Excel Calculations
Easting         Northing        sq m                    sq km
719,510.3358    4,208,764.4633
243,900.3520    4,432,069.0566
--------------------------------------------------------------------
475,609.9837    223,304.5933    106,205,894,001.5610    106,205.8940

Sites used to convert to UTM

http://www.latlong.net/lat-long-utm.html

http://www.rcn.montana.edu/Resources/Converter.aspx

http://www.earthpoint.us/Convert.aspx

Running the same calculation of the difference and the area gives me the result of 106,205 sq km.

The code currently calculates the lat long by doing the following on a polygon of the bbox ...

SqlGeography geog = SqlGeography.STGeomFromWKB(aoi.STAsBinary(), (int)aoi.STSrid);
area = (double)geog.STArea();

This gives me about 9,600 sq km.

What I need to know...

I am trying to find out which calculation is correct.  I am leaning toward 3857 being correct as I the calculations give me approximatly the same numbers no matter how I do them.  Since I feel that 4326 is incorrect, I would need a better way to calculate this in the current application (uses C# and SQLGeometry).  If someone sees an error in my calculations please let me know and how to do this better if so.


texts:
I try to fill the fields with the condition if / else but it does not work.
I do not know where is my mistake in the code.

def TextValue ( !Classe! ):
if !Classe! == "Steppe arbustive"
  or !Classe! == "Steppe arbustive à arborée"
  or !Classe! == "Steppe arborée"
  or !Classe! == "Savane arbustive"
  or !Classe! == "Savane arbustive à arborée"
  or !Classe! == "Savane arborée à boisée"
  or !Classe! == "Savane boisée"
  or !Classe! == "foret galerie relictuelle" :
           return "Végétation naturelle"

elif !Classe! == "Culture pluviale / jachère"
 or  !Classe! == "Culture maraichère" :
           return "Zone de culture"

elif !Classe! == "Cours d'eau"
 or !Classe!== "Mare" :
           return "Surface hydrique"

Codeblock
TextValue ( !Classe! )


texts: Error with CellStatistics arcpy python script

texts:
I have about 20 raster files in a directory. I want to sum all the rasters into one new raster. The code I wrote was:

currentDir=os.getcwd()
tempDir=currentDir+'\\tmp'

arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True

arcpy.env.workspace=currentDir+'\\tmp'
sumRas=arcpy.ListRasters("rr*", "TIF")

outCellStatistics = CellStatistics(sumRas, "SUM", "DATA")

When I run it I get this message:

Traceback (most recent call last):

File "D:\SAPROGtrab\CPIR.py", line 50, in outCellStatistics = CellStatistics(sumRas, "SUM", "DATA")

File "D:\programas\ArcGIS\Desktop10.6\ArcPy\arcpy\sa\Functions.py", line 3173, in CellStatistics ignore_nodata)

File "D:\programas\ArcGIS\Desktop10.6\ArcPy\arcpy\sa\Utils.py", line 53, in swapper result = wrapper(*args, **kwargs)

File "D:\programas\ArcGIS\Desktop10.6\ArcPy\arcpy\sa\Functions.py", line 3169, in Wrapper [function] + Utils.flattenLists(in_rasters_or_constants))

RuntimeError: ERROR 999998: Unexpected Error.

If I run it with only 3 raster files it work ok.

Any idea about this error?


texts: Work order: project then resample, or resample then project

texts:
I'm working with 10m DEMs in northern New York, along the Canadian border (in ArcGIS 10.2). To complete a regional elevation map (and eventually slope, flow, etc.), I need to also include a portion of Canada just over the border. While the USGS provides 10m DEMs for most of Canada's border, there is one small section that only has 30m resolution available.

I'm trying to resample the 30m DEM to 10m resolution, but am not sure which order I should work in. The DEMs are in GCS_NorthAmerica_1983, and I'm projecting them to UTM 18N.

Should I project the 30m DEM to UTM 18N, and then resample? Or should I resample, and then project to UTM 18N?

Any suggestions or thoughts on which order would provide more accurate results?


texts: NO values of field in attribute table is seen in ArcGIS Desktop?

texts:
I am having a problem in my ArcGIS 10 project.

I had download land use data from my city but while opening attribute table the values for one of the fields (Land use values) are not visible.The file was in Shape format

I Checked the fields in the properties, then also no values are visible in that field.

Whenever I turn on the edit mode, the values are visible.

IF I stop editing mode the values disappear.

While classifying the values from symbology for creating a map, the Map is also not visible. I think the map is created but its not visible.
IfI open the file in Q-GIS then all the fields in attribut table can be seen. I feel comfortable working in Arc gis .


texts: Geoserver ( Java.exe ) cause server to run 100%

texts:
Our server runs at 100% CPU when using map services build on geoserver + postgresql
I can not find anything in geoserver logs.

It is Java.exe takes 100% CPU and because of that server can not do
anything else until we restart geoserver.

Environment: Geoserver 2.5.1, PostgresSQL 9.3, SQL Server 2008 R2

Any suggestion what issue causing it to go 100%?

Note/Solution: It was causing when we tried to project huge database native web Mercator into wgs84. we changed the source data back to wgs84 and it never went to 100% again


texts: copy a column from a table to an other using qgis

texts: i have a shapefile and I want to copy two fields from it to an other table wich is empty. If you an idea about how to do it please tell. I'm using qgis 2.18.1.


texts: Make curved street junctions in a cross styled junction in ArcGIS

texts:
I have to make new city blocks based on buildings in a given area. The blocks have to be rectangular and their sides have to be parallel. Then I have to make this automated, so I have to convert that to Python, bad that's not a problem. I just can't fix something by hand.

I used :

2 Minimum Bound Geometry-Envelope, one to group the buildings in the block ( must have a common field), and one with all the building, to take the project area
Converted both to lines
Generated points along lines that group invidual city blocks
Created a Thiessen polygon for these points, and dissolve the polygon based on their common field that is the number of the city block, and got the midline of the cityblocks, which is practically the centerline of road network of the area and turned that into line.

My problem is that in polygons edges the line of the roads makes a curve. Does anyone know how to delete the curved line and make lines straight, with an angle of intersection around 60-90 degrees?
I want to solve that cause I am going to make the city blocks by buffering the road axis to make my cityblocks.
With these curved road segments, my final polygons look like this :

TLDR: How can I make junctions from state 1 to state 2?

I can work on ArcGIS Desktop Advanced and ArcGIS Pro Advance.


texts: Showing Arrow connection relation between Polygon centroids to points in ArcGIS Desktop

texts:
I am looking for a similar solution to Making arrow from polygon centroid to point using QGIS but using ArcGIS Desktop.

I have a polygon and a point shapefile. The polygon area drains to the point node (This is a drainage node). The connectivity is coded in the attribute table name field in shapefile.

Now, I want to visually represent this connectivity using arrows from polygon to point. I also want this arrow indicator to be dynamic.

Meaning, any change in the attribute table reference will update the arrow connectivity. All these features are already available in QGIS. I want to replicate the same mechanism in ArcGIS Desktop.


texts: How can I use pyqgis to only show point that do not share a common value?

texts: I have two layers of point data containing the coordinates of addresses from two different years. I want to display only the new addresses on the most recent dataset that do not exist in the database from the earlier year. To do this I am trying to use a feature which both layers contain, the value of which is a unique identifier for each address which is the same across years. I want to know how I can use pyqgis to display only these new addresses by making a script which will look through the feature in question on both datasets, and if the value of said feature on one dataset is the same as the other dataset, neither point the will be displayed. Any suggestions?


texts: SAGA hydrological analysis in QGIS 2.8 on Mac

texts:
I am running QGIS 2.8 on a Mac machine and I want to calculate the catchment area using SAGA - Catchment Area. It starts running but after a while when QGIS wants to open the generated files there is an error message that QGIS cannot open the files:

Oooops! The following output layers could not be open

Name for output filled (flooded) elevation raster map:
/var/folders/g1/r0fk8c2j5ddf7p5jf_1x851w0000gn/T/processing/b30a24aad2ba451abc6875ff7ec46740/filled.tif

Name for output flow direction raster map:
/var/folders/g1/r0fk8c2j5ddf7p5jf_1x851w0000gn/T/processing/93e5878769d84ebc981bfb3721b1a861/direction.tif

Name for output sink-watershed raster map:

In contrast the  SAGA - Channel Network and drainage basin analysis fully works. What is wrong here?


texts:
As written in the Title: How do I group layers in QGIS?
I have tried dragging them on top of each other and looked in the layers setting but could not find anything.


texts: Does SDE.ST_LENGTH() calculate the true length of arcs?

texts:
Oracle 18c/10.7.1 EGDB
I have a polyline FC that has features with arcs.
Does the Oracle SDE.ST_LENGTH() function calculate the true length of arcs?
Or does it only calculate the straight-line length between the arc's vertices?


texts:
I am new to QGIS. My goal has been to convert a .shp file to a raster, eventually to a .bmp for an input map into a separate piece of software I am using... this separate software also needs the map to be in UTM.
My steps after inputting vector layer:

Change CRS from Conus Albus to WGS 84 / UTM.
Use Rasterize... Conversion- Vector to Raster (using Byte as the output data type)
Raster...Conversion- Translate (Convert Format)- Save to File- BMP

When I open this .bmp from my desktop (viewing in XnView or Windows Photos) the .bmp is severely tilted to the right...
I must be missing something simple? Any advice?

Second photo is the output.


texts:
When trying to save a style to postgis database through QGIS I am getting the error:

  'Unable to save layer style. It's not possible to insert a new record into the style table. Maybe this is due to table permissions (user=admin_gis). Please contact your database administrator.

After reading other posts, I have found there is no layer_styles in the public schema. I have tried adding this table in manually and re-saving but still receiving the same error.

I cant seem to find a resolution that works for me online, has anyone else had/solved this problem?


texts:
Is there a difference between the actual zooming when using the setView and panTo methods in Leaflet?

map.setView([52, 13]);

map.panTo([52, 13]);

The zooming effect looks identical to me, or am I missing something? Are they the same unless you use the optional parameters?


texts:
In part of a city, I am trying to count rooftops on a basemap/Google satellite in QGIS and create several polygons around x number of rooftops in the area. One similar question like mine was previously asked, but the response was only explained how to do it in Arcmap, not QGIS.

Is there anyway to do it in QGIS?


texts:
I have multiple WMS layers. One of the layers is from a local GeoServer (UAV generated GeoTIFF). I am unable to zoom GeoTIFF up-to it's full resolution.

How to enable deep zoom in to WMS above 20% zoom level?

Here is my code:

<!DOCTYPE html>

    <title>Sample Web GIS using Leaflet and Geoserver</title>
    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1.0">
    <link rel="shortcut icon" type="image/x-icon" href="docs/images/favicon.ico" />
    <link rel="stylesheet" href="D:/Web_GIS/html/leaflet/leaflet.css"/>


texts: ArcGIS create profile graph from existing point feature

texts:
I have a dataset of X,Y points (with a spacing of 20m between the points) that represents a flowline along a glacier. Next to dat, I got another point dataset of the glacier height with X,Y points and Z data for the elevation. What I want to do is extract the surface profile of the glacier along the flowline. Is there a way this can be done?

Attachment:
brown points = flowline (X,Y)
blue points = glacier shape with elevation (X,Y,Z)
TIN surface of the glacier elevation at the background


texts:
I am using this:

COPY data_table (col1,col2,col3) FROM 'datafile.csv' DELIMITERS ',' CSV HEADER;

to get csv data into my postGIS table. It works wonderfully, except I need to have a 'pk_id serial PRIMARY KEY' in the database table as well, which will then make the columns not match up between 'data_table' and 'datafile.csv'. And I think this will give an error when I run the COPY command.

I could add the pk_id serial PRIMARY KEY after the import, but I want to run many imports into data_table so this doesn't seem like it would work.

So my question is: what is an alternative way to get csv data in - can I still use COPY? Is there a way like in mySQL to "INSERT INTO table_name (col1, col2,col3) VALUES (val1,val2,val3) where I specify the affected database columns?

Thanks


texts: Create a closed polygon from a series of edges given a point

texts:
I have a series of points fed into a voronoi graph class that returns a series of edges defining the bounds between these points, the process works it returns a series of edges.

I now need to do a ray/fill using the points and these edges to build polygons around each point.... for the life of me I can't seem to figure out a process of doing this or a code sample...
C# code or thought process would be great .... maybe ray collision?


texts: Find Largest Building in each Parcel (largest feature in each groupby set)

texts:
I have a question similar to this one, but the proposed approach doesn't seem to work.
For context, I am new to qgis.
Given a set of parcels, each with 0 or more buildings, I want to select only the largest building in each parcel.
There are about 200k parcels and 320k thousand buildings.
Example
Building A:
 parcel_APN: 1 (ParcelNumber)
 SHAPE_Area: 100

Building B:
 parcel_APN: 1
 SHAPE_Area: 500

Building C:
 parcel_APN: 2
 SHAPE_Area: 50

>>> LargestBuildings: [B, C]

I tried two approaches:
SQL Query DB Manager
Seems like this should work, but I let it run for 2 hours and it did not finish.
Dataset is big, but seems to be shouldn't take that long for a SQL join.
A simpler SELECT with MAX area and groupby but without the join take less than  5 seconds.
The approach is to first group by APN to find the largest shape area, then do inner join using id and max area , with the same table, to retrieve the largest buildings.
SELECT * FROM
 (
   SELECT
     parcel_APN,
     MAX(SHAPE_Area) as max_shape
   FROM buildings as t0
   GROUP BY parcel_APN
 ) as t1
JOIN buildings as t2
ON t1.parcel_APN = t2.parcel_APN and max_shape = t2.SHAPE_area

Using Attribute Expression
This is the approach suggested in the SO mentioned above:
maximum("SHAPE_Area", group_by:="parcel_APN")

I let this run last night, this morning was still running and eventually QGIS crashed.
Questions
Do these approach seem correct? Is there an easier way to go about this?
Are the file sizes too big for something like this in QGIS? And maybe I need to work in subsets or work with the data outside?

Update
I created a small subset to test, and the attribute expression method does not seem to work - selects all buildings:

Update 2
DB Manager Query seems to work! I guess now I just need to find a way to run it without crashing QGIS, or load it into PostGIS to run it.


texts: Unable to process batch buffer using ArcPy

texts:
I'm trying to buffer multiples files at once using ArcGIS. My code looks correct until the last line where it returns error related to buffer tool. I've used this code.
import arcpy
import os
from arcpy import env
path = "G:\\Python\\Arcgis\\Materials"
env.workspace = "C:\\Users\\XYZ\\Documents\\ArcGIS\\Default.gdb"
arcpy.env.overwriteOutput = True
inputlist = []
shapelist = os.listdir(path)
shapelist = [x for x in shapelist if x.endswith(".shp")]
for shape in shapelist:
    print(shape)
    inputlist.append(shape)
output = "G:\\Python\\Arcgis\\Buffer"
distanceField = "100 Feet"
arcpy.Buffer_analysis(inputlist, output, distanceField)

Also, I've tried it in ArcGIS Desktop as well as ArcGIS Pro. The files are projected.
The error while using ArcGIS Desktop is as follows:
PS G:\Python\Arcgis> & c:/Python27/ArcGIS10.4/python.exe g:/Python/Arcgis/buffer_standalone.py
area1.shp
area2.shp
area3.shp
Traceback (most recent call last):
  File "g:/Python/Arcgis/buffer_standalone.py", line 15, in <module>
    arcpy.Buffer_analysis(inputlist, output, distanceField)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\analysis.py", line 692, in Buffer
    raise e
RuntimeError: Object: Error in executing tool


texts: Defining projection surface height in proj4 string

texts:
I'm using proj4 in an actual project, which needs accurate projection calculation. To realize accurate projection, it defines a projection surface height, which is calculated from the average altitude of the working zone. It is like, for example:
If the projection surface height is 100m, and I'm using WGS84 datum, I have to add 100m to an of the datum and maintain the eccentricity(e) of the datum, which will make the b become (a+100)(1-f), where f=1/298.257222101
Now my working solution is directly using the a and b computed as mentioned above in the string, like
+proj=tmerc +lat_0=0 +lon_0=30 +k=1 +x_0=500000 +y_0=0 +a=6378237.0 +b=6356852 +units=m +no_defs

where +a=6378237.0 means the a of WGS84(6378137.0) plus 100m, +b=6356852 means 6378237.0*(1-1/298.257222101), f is the oblateness of WGS84 datum.
Now I'm wondering is there any parameter in the proj4 string that can let me directly specify the projection surface height directly, instead of manually giving a and b? like +proj_surface_height=100m?


texts: Largest polygon inside a GeoJSON multipolygon

texts:
How to find the largest polygon in GeoJSON of type multipolygon so that I can add a Leaflet tooltip to the largest polygon rather then a random polygon in multipolygon?

I have a Multipolygon and when adding a tootltip to it with Leaflet, the tooltip is attached to the first polygon or a random polygon in the multipolygon. I don't want it to be that way. I want to add tooltip to the largest polygon.


texts: Defining grid-names: Setting up a 1m2 grid over a 100x100 meter area

texts:
I'm working on an old archaeological field-survey dataset from the 1960's. I'm trying to create a polygon grid (fishnet) over a specific area where I can give the individual 'tiles'of the grid an address as attribute (for example grid square Aaa-1). I need this 'address' in order to link the specific squares with  the data of the finds attached to this square (determinations, types, quantities, etc).

About the definitions of 'addresses': Finds are recorded in 1m2 squares within a grid of 100x100m meter. Along the x-axis the squares are numbered numerical 1-100 and along the y-axis alphanumerical (for example: Aaa-Aaj for the first 10 rows, Aba-Abj for the second 10 rows, etc until Ajj, which is nr 100.).

I've  tried the FSC tools, with the OSGR Tool to create a grid, but this didn't worked out due to an inclination in my grid and it is not possible to position the grid on specific coordinates. The 1966 grid was not orientated to the true north, but set-out with an angle of 5,5 degrees in regards to the north. I have the coordinates of the four corners of the 100x100 meter square.

I'm currently working in QGIS, but if (easy) solutions are possible in other GIS environments, I'd be happy to hear.


texts: R: Overlay using functions containing if/else statements

texts:
I have two rasters, x and y, and I'd like to do calculations with them using the overlay function. For example,

MyFun <- function(x, y) x * y
res <- overlay(x, y, fun = MyFun)

This seems to work fine. Now, I want to add a conditional statement to my function like this:

MyFun <- function(x, y) {
    if(x > 0){
        x * y
    }else{
        abs(x) * y
    }
}

This gives a warning message like below:

Warning messages:
1: In if (x > 0) { :
  the condition has length > 1 and only the first element will be used

Clearly, R isn't happy with my code.

How should I write conditional statements in my functions so that they are compatible with overlay and rasters?

NB This is a toy example. My real task is more complicated, but the principle is the same.


texts: Compute Shortest Distance(s) between two points

texts:
I am using Mapinfo 11. I want user to be able to find shortest route between two points so that user can travel from point A to point B using these path. Any Idea how to achieve this?

Edit: I need to be able to do this with the Bing Hybrid Map.


texts: How to take cosine squared of raster using python

texts:
I have an equation y = cos^2(x) where x is a raster.

How can I write this equation in python?


texts:
I follow the instructions here very closely and yet every time QGIS gets caught in a loop. It takes a few small rasters totally less than 100mb and creates HUGE gb files and will not stop until python.exe is killed in task manager.

From one of the ecw's metadata :

Driver
GDAL provider
ECW
ERDAS Compressed Wavelets (SDK 5.0)
Dataset Description
E:/work/1680_U/aerials/1680_naip/global_mapper/defbri_gm/500ft_clips/87.ecw
COLORSPACE=MULTIBAND
COMPRESSION_RATE_TARGET=1
VERSION=2
Band 1
Band 2
Band 3
Band 4
Dimensions
X: 36702 Y: 14147 Bands: 4
X : 18351,Y 7073
X : 9175,Y 3536
X : 4587,Y 1768
X : 2293,Y 884
X : 1146,Y 442
X : 573,Y 221
Origin
614731,4.61751e+06
Pixel Size
1,-1
No Data Value
*NoDataValue not set*
Data Type
Byte - Eight bit unsigned integer
Pyramid overviews
Layer Spatial Reference System
+proj=utm +zone=16 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Layer Extent (layer original source projection)
614730.9328588016796857,4603364.5349677717313170 : 651432.9328588016796857,4617511.5349677717313170
Band
Band 1
Band No
1
No Stats
No stats collected yet
Band
Band 2
Band No
2
No Stats
No stats collected yet
Band
Band 3
Band No
3
No Stats
No stats collected yet
Band
Band 4
Band No
4
No Stats
No stats collected yet


texts: Clip_analysis producing an extra output in Arcpy?

texts:
I am trying to clip a feature class of floodplain polygons using another polygon object in arcpy, like so:

arcpy.Clip_analysis(in_features = floodplain_polys,
                    clip_features = <PolygonObject>,
                    out_feature_class = path.join(scratch_gdb, 'test'))

The clip runs fine, but it creates two outputs: one called test, and another called floodplain_polys_Clip, both in the same fgdb. I know that the second name is the default output name for a Clip analysis, but I'm struggling to figure out why it's created in the first place. Is this happening because I used a Polygon Object instead of another feature class? If so, why?


texts: Deleting vector lines beyond image extents using QGIS?

texts:
I am using QGIS 2.2

I have one vector line layer (roads) that stretches beyond my base image (see below).

I wish to have the roads "clipped/trimmed or deleted" right at the edge of the image's extent (see below):

How can I achieve this?