texts: KMZ file in wrong coordinate system and won't transfer into correct one

texts:
I'm using a pyroclastic flow simulation model (Titan2D Mass Flow Simulation GUI) that uses UTM WGS 1984 and from what I can tell, the coordinates that I'm inputting are all correct. When I view the output files within the program, they're completely correct (starting points, extents etc all use coordinates and all are correct) but when I put the downloaded KMZ files into Google Earth or any GIS program, they end up in UTM Zone 31N instead of in Tenerife, where I want them. They're at the correct latitude but are just in Algeria instead of in Tenerife!
I have access to QGIS, ArcMap and ArcGIS Pro and have been using ArcGIS Pro to try to 'Define Projection' and 'Project' and anything else I can think of (turned them into shapefiles and tried again to 'define projection' and 'project') but nothing seems to be working at all.
Is there any way to move them along the line of latitude from Algeria to Tenerife?
Does anyone have any suggestions as I need to move them over to Tenerife so I can create hazard maps for my dissertation?


texts: Executing buffer analysis in ArcPy with for loop?

texts:
I'm attempting to run the Buffer Analysis tool with a distance of 3 km over multiple shapefiles at the same time. I was given a folder with different types of files and ran a for loop to create a new list with just the files ending with ".shp" that I have below.

import os
import arcpy
from arcpy import env
env.workspace = "file_path"
env.overwriteOutput = True
items = os.listdir("file_path")
 newlist = []

for names in items:
    if names.endswith(".shp"):
        newlist.append(final)
print newlist

Thus, this creates my list with just the shapefiles. However, this is where I struggle. I'm not sure how to run the buffer analysis tool over this new list that I've created.


texts: Visually show input data from attribute table on map

texts:
I have 8500+ point-locations in QGIS, all with unique ID numbers, showing on my map; however, the unique ID numbers are not displayed.

Is there any way to filter this so that the ID numbers are present next to the point (all points at the same time)?

See pictures for clarifications.

https://i.ibb.co/3kLqnvY/avs-id.png


texts: How to write a single pixel or a line of pixels to GeoTIFF using GDAL/Python

texts:
Let's say I want to write a single pixel or a 1D array of pixels (e.g.all at the same latitude) of pixels to a single-band of a Float32 GeoTIFF using the Python GDAL bindings.  How can I do so?  Or is it not possible?

I currently have the following method for truly 2D ndarrays. But occasionally I get requests for a pixel or line of pixels, which cause this method to fail.

raster_band.WriteArray(my_pixel_array)

FURTHER CLARIFICATION:

I'm trying to write a new image where the entire extent of the new image is either just one pixel or just a line of pixels. Thus, I've tried just setting xoff and yoff to 0, and this fails the same way as without xoff and yoff. The error I'm getting is expected array of dim 2, e.g. when the shape of my ndarray is (4,) (and when the GDAL raster is initialized with RasterXSize to 4 and RasterYSize to 1 via a pertinent out_raster.SetGeoTransform (EPSG:4326 location and size) method).


texts: Viewing Tables in the AxTOControl using ArcObjects

texts:
I was wondering if anyone out there has written a Table filter for the AxTOCControl in the ArcObjects Library.

I've rolled my own Table view for now in C#, but would like to try and keep things consistent.

If you could point me to any resources I would be very appreciative.


texts:
I want to use present year in a CLASS in a map file. I managed to write permanent expressions like:

LAYER
    CLASSITEM "YEAR"
...
    CLASS
        NAME "Previous years"
        EXPRESSION ('[YEAR]'<'2015')

But I  failed to construct a working expression, which includes a reference to present time.

There is an example of what I would like to write:

LAYER
CLASSITEM "YEAR"
...
    CLASS
        NAME "Present year" # Year n
        EXPRESSION ('[YEAR]'=='([date format=YYYY]'))
...
    CLASS
        NAME "Previous year" # Year n-1
        EXPRESSION ('[YEAR]'=='([date format=YYYY]'-'1'))
...
    CLASS
        NAME "Archived years" # Years <n-1
        EXPRESSION ('[YEAR]'<'([date format=YYYY]'-'1'))

I adapted [date format=YYYY] expression from a MapServer template page.It didn't work.... But unfortunately I did find no simple examples of templates imported variables direct use in a map file, especially for [date] variables.

Has someone any idea of what it could look like?
I'm just looking for reusable examples that I could adapt for my scripts.


texts: Topographic Wetness Index, Raster Calculator Error

texts:
I'm having issue calculating TWI, my projected coordinate system is NAD_1983_UTM_ZONE_17N. As I'm in Ohio, USA so I'm using this coordinate system. I have calculated the FA, FD and Slope but when ever I try to calculate TWI from these values using raster calculator, it shows an error. My other values are:
Columns and Rows: 1408, 1188
Number of Bands: 1
Cell Size (X, Y): 1, 1
Uncompressed Size: 6.38 MB
Format: FGDBR
Source Type: Generic
Pixel type: floating point
Pixel Depth: 32 bit
Pyramids: Level:3, resampling: Nearest Neighbor

Note: Attached is the images also.


texts: Creating TIN, Surface Difference python script with ArcPy?

texts:
I have been trying to develop a master script to go through a series of geoprocessing tools and executing them. I am trying to generate a TIN and contours from a raw xyz file.  Below is a very simple script that integrates the various tools. I am having difficulty inserting the correct inputs and outputs for the TIN generation and Surface Difference calculation/generation. Let me know if i am missing an input/output specification or a simple default value. I tried doing this same workflow on Model Builder to then export its python script, however i was unable to execute a successful TIN Output since for some reason only allowed for a Raster Output.

Note: I decided not to insert variables, i went directly to its filename.

import arcpy, exceptions, sys, traceback
from arcpy import env
from datetime import datetime

start_time = datetime.now()
print "The program has began processing", start_time

try:
    #Enables ArcGIS 3D Toolbox
    arcpy.CheckOutExtension("3D")
    #Set workspace
    env.workspace = "C:\Standalone"
    #Set Spatial Reference using Well-Known ID  NAD 27 State Plane Missouri East FIPS
    sr = arcpy.SpatialReference(26796)
    #XYZ Processing to Point Shapefile
    arcpy.ddd.ASCII3DToFeatureClass("004d1214_3x3.xyz", "XYZ","004d1214_el.shp", "POINT", "1",sr,)
    arcpy.AddMessage("Processing Complete")

    #Add XY Coordinates
    in_features = "004d1214_el.shp"
    arcpy.AddXY_management(in_features)

    #Create TIN surface out of the generated point shapefile
    arcpy.ddd.CreateTin_3d("tin004d1214_el",sr,"C:\Standalone\004d1214_el.shp POINT_Z Mass_Points <None>","DELAUNAY")

    #Delineate TIN Area. Verify that maximum edge length is greater than the point spacing.
    arcpy.ddd.DelineateTinDataArea("tin004d1214_el", 10, "ALL")
    arcpy.AddMessage("Complete")

    #Trigger Surface Difference to yield new elevation, Insert optional directories for TIN output
    inSurface = "tin004d1214_el"
    inReference = "CLASSIFIED PATH DIRECTORY"
    outTin = "tin004d1214_lrp"
    outPoly = arcpy.CreateUniqueName("004d1214_lrp.shp")
    arcpy.SurfaceDifference_3d(inSurface,inReference,"C:\Standalone\004d1214_lrp.shp",outTin) #Check for syntax

    #Generate Surface Contours off the Reference Survey surface
    arcpy.SurfaceContour_3d(outTin,"004d1214_contours.shp")  #Check for contour field add

except arcpy.ExecuteError:
    print arcpy.GetMessages()
except:
    # Get the traceback object
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    # Concatenate error information into message string
    pymsg = 'PYTHON ERRORS:\nTraceback info:\n{0}\nError Info:\n{1}'\
            .format(tbinfo, str(sys.exc_info()[1]))
    msgs = 'ArcPy ERRORS:\n {0}\n'.format(arcpy.GetMessages(2))
    # Return python error messages for script tool or Python Window
    arcpy.AddError(pymsg)
    arcpy.AddError(msgs)

end_time = datetime.now()
print ('Duration: {}'.format(end_time - start_time))


texts: Openlayers 3. OnClick on touch device not triggering perfectly

texts:
I have the forEachFeatureAtPixel function and try to trigger an event on click on a point layer with a custom icon. However, the performance on touch devices is really bad. I need to click several times until it triggers the event (sometimes it doesn'tdo it at all..) . I also used the moveTolerence feature for the map, but still have the same issues. On mouse click on my desktop PC everything works like a charme. What am I missing to make the performance as smooth as on a desktop PC?

This is my code:

$(document).ready(function() {
map.on('singleclick', function(evt) {
    var feature = map.forEachFeatureAtPixel(evt.pixel, function(feature, layer) {
        return feature;
    }, {hitTolerance:10});
    var layer = map.forEachFeatureAtPixel(evt.pixel, function(feature, layer) {
        return layer;
    }, {hitTolerance: 15});
    if (feature) {
        if (layer === vectorLayerFB) {
            window.open('https://www.example.com', '_blank');
        }
        if (layer === vectorLayerBlog) {
            window.location = 'www.example.com';
        }
    }
});
[...]
});


texts:
I am working in c++ qt4.8.5 and api qgis 2.8.2 c++.
In my application I had to put polygon to define attractive zone.
I found a great post in stackexchange (86812).
I have based my on this post and on the pyqgis cookbook,
but in c++ it seems that it is more difficult.
I create a layer in memory with a polygon type,
then i create a function to generate ( for the moment ) a random polygon shape with 10 points.
There is no error in the build, when i call the function there is nothing on the display.
I check the coordinate of the point it is ok.

here is my code

create the polygon layer :

QgsVectorLayer * lVectorLayer = new QgsVectorLayer("Polygon",pType,"memory");
lVectorLayer->isValid();
qDebug() << "Vector valid :" << lVectorLayer->isValid();
// add to the map
QgsMapLayerRegistry::instance()->addMapLayer(lVectorLayer);
// insert the vector layer into the vectorlayermap
mMapLayer.insert(pType, lVectorLayer);
// Add the Layer to the Layer Set
QgsMapCanvasLayer PolygonLayer = QgsMapCanvasLayer(lVectorLayer, TRUE);
mpS57->myListForPrint.insert(1,PolygonLayer);
// Set the Map Canvas Layer Set
mpMapCanvas->setLayerSet(mpS57->myListForPrint);

then the function to create the random polygon :

QgsVectorLayer * lLayer = mMapLayer.value(pType);
if (lLayer == NULL) qDebug() << " lLayer NULL";
else qDebug() << "lLayer Not NULL";

qDebug() << "Vector is Valid ? ::" <<   lLayer->isValid();

QgsVectorDataProvider * lDataProvider = lLayer->dataProvider();

QgsFeature MyPoly = QgsFeature();

QgsPoint lPoint2;
QgsPolyline lpolyline;
QgsPolygon lPolygon;
// origin point
QgsPoint lOrgPoint;
QgsPoint lPoint1;
PosX =(qrand()%(Xmax-Xmin)+Xmin);
PosY =(qrand()%(Ymin-Ymax)+Ymax);
PosX=PosX / 100000;
PosY=PosY/10000;
lOrgPoint = QgsPoint(PosX,PosY);
lPoint1 = lOrgPoint;
// add 10 points
for (int i=0; i<10; i++) {
    //calc new random X and Y
    PosX =(qrand()%(Xmax-Xmin)+Xmin);
    PosY =(qrand()%(Ymin-Ymax)+Ymax);
    PosX=PosX / 100000;
    PosY=PosY/10000;
    lPoint2=QgsPoint(PosX,PosY);
    // create the line
    lpolyline.append(lPoint1);
    lpolyline.append(lPoint2);
    // insert the line in polygon
    lPolygon.append(lpolyline);
    // clear the line
    lpolyline.clear();
    //replace lpoint1 by lpoint2
    lPoint1 = lPoint2;
}
// close the polygon
// create the line
lpolyline.append(lPoint1);
lpolyline.append(lOrgPoint);
// insert the line in polygon
lPolygon.append(lpolyline);
// clear the line
lpolyline.clear();

MyPoly.setGeometry(QgsGeometry::fromPolygon(lPolygon));

QgsFeatureList lPolyList;
lPolyList << MyPoly;
lDataProvider->addFeatures(lPolyList);
lLayer->updateExtents();
QgsMapLayerRegistry::instance()->addMapLayer(lLayer);

as you can see, in c++ is more complex than in python
i try to put a list of point in the frompolygon(listofpoint) but there is an building error in c++ we had to respect the type, in this a QgsPolyline.

I don't understand why tehre is nothing on the display
there is the update, the layer is added to the map, the feature is added to the dataprovider, everything seems to be right but nothing.

any idea ?

added information 09-10 : I have made other test with a simple line, same result, but with the same list of point if I use the qgsrubberband, the rubberband is displayed fine on the map, the test with the rubberband confirm that the list of point are correct.


texts: Selecting groups of files and iterating composite raster

texts:
I have a query regarding the best way to process the multispectral data I requested to be collected for my MSc study site.

The data I am working with has been saved as single band tiffs, rather than as multiband ones. They are named as follows.

1454327096.tif = Green band
1454327096_1.tif = Red band
1454327096_2.tif = Red_edge band
1454327096_3.tif = NIR band

I have about 1084 of these images so I would like to find some way of automating the process of compositing the bands.

My thought was to use an iterator in modelbuilder in ArcGIS. However I have no idea how to get the iterator to group files together and only composite the ones where the file names are the same.


texts: Where I should add cql_filter param in WMTS layers by OpenLayers v6?

texts:
I am going to use cql_filter in wmts layers,this is my code, I read OpenLayers v6 Document , but I couldn't find params, I googled but I could'nt find any soltion, about cql_filter in WMTS Layers, If I want use WMS Layers I should insert "cql_filter" in ol.source in params option, but WMTS layer has not params, Hence, Where I should insert "cql_filter" If I use params in my code I add params like this
  params: {'LAYERS': `${geoServerWs}:${layers_}`,'CQL_FILTER':"objectid=16"},

to my code but it didn't show feature based on cql_filter
    var wmtsLayers=[
        new ol.layer.Tile({
          title:`${layers_}`,
          opacity: 1,
          source: new ol.source.WMTS({
            //attributions: '...',
            url: `http://localhost:8080/geoserver/gwc/service/wmts?`,
            layer: `${geoServerWs}:${layers_}`,
  params: {'LAYERS': `${geoServerWs}:${layers_}`,'CQL_FILTER':"objectid=16"},
            matrixSet: 'EPSG:900913',
            format: 'image/png',
            serverType:'geoserver',
            tileLoadFunction:function(tile, src) {
                var xhr = new XMLHttpRequest();
                xhr.responseType = 'blob';
                xhr.open('GET', src);
                xhr.setRequestHeader('Authorization', 'Basic ' + window.btoa(`${geoUserName}` + ':' + `${geoPass}`));
             xhr.setRequestHeader('Content-Type', 'application/json');

                xhr.onload = function() {
                    var objectURL = URL.createObjectURL(xhr.response);
                    tile.getImage().onload = function() {
                        URL.revokeObjectURL(objectURL);

        };
        tile.getImage().src = objectURL;
      };
      xhr.send();
    },
           projection: projection,
            tileGrid: new ol.tilegrid.WMTS({
              origin:ol.extent.getTopLeft(projectionExtent),
              resolutions: resolutions,
              matrixIds: matrixIds,
              extent: ol.proj.transformExtent([44.03, 23.34, 63.34 ,39.78], 'EPSG:4326', 'EPSG:900913')

            }),

          })
        })
      ];


texts: What Professional / Industry Associations are available that are related to GIS?

texts:
This was a well liked (but not hugely rated) question in the beta, so here goes.

Please provide a brief synopsis and links if possible. I'd also like to hear about things that are in the planning stage.

What group or association do you belong to? What are you working towards? Why?


texts:
Does anyone know any good freeware alternatives to eCognition?

I need to do image segmentation and classification. I tryed the ecognition trial and it was really good. I'm doing the job with SPRING, but it has a lot of limitations...

Is there any qgis plugin?


texts: I have a set of several tiff:s, the largest is 4,7 MB, that I want to create a TMS out of. If I don't want to use Gdal and write commandlines is there a free software option that I can use?


texts: Given two Shapefiles, how to tell which one's CRS is set correctly?

texts:
I have a polygon shapefile projected in greekgrid (GGRS87_Greek_Grid)
then I open an other shapefile which is in ETRS_1989_LAEA
In ARCGIS 9.3 the two shapefiles are displayed correctly.
But if I do the same in QGIS there is a shift of the ETRS_1989_LAEA shapefile

How can I find which of them is correct?


texts:
I've implemented labelling in my code and solved a fair amount of problems that seemed insurmountable at first. Now that the labelling is working I have directed my attention to the positioning of the labels and am extremely frustrated with the capabilities offered by the base ArcEngine.

I have no experience with Maplex (or any other tool) and suspect that that extension (or a 3rd party tool) is the only way to do decent labelling in the ArcGIS product line. Is this a true picture of the general experience? If so, what solution would you, dear reader, offer? Is Maplex the best and only way? Are there acceptable 3rd party solutions?


texts: LAStools blast2dem: error writing data for field "StripOffsets"

texts:
I was processing an overage file to generate DTM and DSM using blast2dem. The generated DTM file can be viewed using ArcMap and ENVI. However, ArcMap and ENVI cannot open the generated DSM file. I don't know exactly what is the cause of this problem but there is an error reported every time I try to repeat the same process of generating DSM. It says:

  Error writing data for field "StripOffsets".

Please see the snapshot below. Does anyone know how to resolve this issue?


texts: OpenLayers 4 Displaying GeoJSON Point Data

texts:
I'm struggling to properly display my GeoJSON data which is being serialized from a GeoDjango model, though hopefully that isn't relevant to the issue at hand.

I copy and pasted a GeoJSON example from the OL4 docs and then tried to add in my own GeoJSON just to try to get it displaying properly but I have so far been unsuccessful.

Here's the broken jsfiddle to show my problem. Where am I going wrong here?


texts: Why there is no "View Output" button for the command "v.to.points" in Grass plugin in QGIS?

texts: How do i get or where do i get my output vector point layer after i run the command v.to.points to my vector line layer?


texts:
I have a strange behaviour in QGIS with GRASS tools. I am trying to create basins in the alp region and use he SRTM3 data. As I can see on the flow direction raster the SRTM3 data seems to act nice:

I was using the r.watershed algorithm with a minimum size of 1000, 2500 or 5000 and I am always seeing the same problem:
straight lines across a slope and not following the main flow direction raster (visible in the center of the above screenshot).

Can anyone give me some advice on how to take care of this issue?

I was using QGIS 2.12.3, Ubuntu OS, the SRTM3 39_03 dataset and lots of coffee.
I tried D8 and non D8 algorithm as well as different values for RAM usage...
I will try it once more with an older version of QGIS, but it seems crazy!


texts: Automatic point plotting in ArcGIS from thresholded image

texts:
I would like to know whether it is possible to plot points automatically on a map in ArcGIS by using a raster file as a reference.

The image below would represent an example of line that comes from thresholding procedures using OpenCV. Each point would be plotted following the white line (in contrast with the black background):

I would have done it manually but the number of lines if huge.

My view would have been to ask ArcGIS to recognize the white pixels positions and affect a coordinate for each pixel. The output would be a table that can be exported in a .csv file.

If someone has a better idea or you think it is not possible, I would be happy to hear your thoughts.

Edit

Only points actually. I though about coding something similar the following.

We extract white pixels only:

for i in range(img.shape[0]):
    for j in range(img.shape[1]):
        if img[i][j] == "255":
            ... to be defined

Then we call the Raster To Point tool from ArcGIS to convert white pixels only into points:

arcpy.RasterToPoint_conversion("image_with_white_pixels.jpg", "c:/output/source.shp", "VALUE")

Then I guess I would define the coordinate of the points using the Add XY Coordinates tool (from Data management toolbox).


texts: "Geographic coordinate systems warning" when add a file in ArcGIS

texts:
I have two files to display in ArcGIS.

First, I add file one, whose Geographic Coordinate System is "GCS_North_American_1983". Then the second file is added, whose Geographic Coordinate System is "GCS_WGS_1984".

At this moment, a Geographic coordinate systems warning occured as shown below.

I know that it's because these two files have different Geographic Coordinate System. My question is :

Does this matter (As it is a "warning", sometimes it doesn't matter)?
If it does matter, what should I do?


texts: QGIS - Set display order of polygon features in the same layer

texts:
in QGIS I have a polygon layer with overlapping features. While this is not ideal (and I am currently changing this), I am wondering whether I can in any way affect the order of the features displayed. My smaller feature is in the background, so I can´t see it (and have a harder time editing it, though I am solving this by making the overlapping layer transparent).

I have tried sorting the attribute table differently, this does not affect the order of display of the features in my layer.

Any ideas?


texts: Creating DEM from contours using ArcGIS Desktop?

texts:
I have a large set of contours that I received broken down into several different subsets. I have merged most of the subsets together creating a larger set of contours in three different files. I have then created three different raster datasets from the contour shapefiles, however in that process, the rasters have a lot of zero elevation areas that shouldnt be there.

How do I get rid of these zero elevation areas to make a more accurate digital elevation model?

I was hoping the topo to raster tool would interpolate in between the different contours to clean up those areas but it did not. But, some of the zero elevation areas are in places where there is a stream, pond, or reservoir.


texts: Road distance travelled between multiple postcodes

texts: Is there anyway of using QGIS offline... not using API... I want to calculate road distance travelled between multiple postcodes to single destination.. it works perfectly for small number (200) but software gets crashed or through Python error as I have 3000+ postcodes to measure... I created shape file using MMQGIS->Geocode and for distance travelled calculation used OSM tool.


texts:
I am using Coursened Exact Matching to improve balance in my sample. The summary matching results after applying CEM are that 11,128 pixels have been matched:
Sample sizes:
          Control Treated
All         19741    8283
Matched      6728    4400
Unmatched   13013    3883
Discarded       0       0

Firstly, is it possible to specify 1:1 matching with CEM? I want to only retain matched pairs, which requires the same amount of pixels in my treatment & control groups.
I am also having trouble getting a unique identifier for each matched pair so that I can calculate the difference between the treatment and control. I am using the code:
pair: Produces a paired sample out of a CEM match solution

When I execute the below command, it says there are 0 matches.
 > # we want a set of paired units

> pairsample <- pair(match.out_cem, data)

Total number of units paired in CEM strata: 0
Total number of units matched: 0

Can anyone help or point me in the right direction?
Thanks, Tori


texts: How do I correctly reproject a geodataframe with multiple geometry columns?

texts:
In the geopandas documentation it says that

A GeoDataFrame may also contain other columns with geometrical (shapely) objects, but only one column can be the active geometry at a time. To change which column is the active geometry column, use the set_geometry method.

I'm wondering how to use such a GeoDataFrame if the goal is to flexibly reproject the geometrical data in these various columns to one or more other coordinate reference systems. Here's what I tried.
First try
import geopandas as gpd
from shapely.geometry import Point

crs_lonlat = 'epsg:4326'        #geometries entered in this crs (lon, lat in degrees)
crs_new = 'epsg:3395'           #geometries needed in (among others) this crs
gdf = gpd.GeoDataFrame(crs=crs_lonlat)
gdf['geom1'] = [Point(9,53), Point(9,54)]
gdf['geom2'] = [Point(8,63), Point(8,64)]

#Working: setting geometry and reprojecting for first time.
gdf = gdf.set_geometry('geom1')
gdf = gdf.to_crs(crs_new)   #geom1 is reprojected to crs_new, geom2 still in crs_lonlat
gdf
Out:
                             geom1          geom2
0  POINT (1001875.417 6948849.385)  POINT (8 63)
1  POINT (1001875.417 7135562.568)  POINT (8 64)

gdf.crs
Out: 'epsg:3395'

So far, so good. Things go off the rails if I want to set geom2 as the geometry column, and reproject that one as well:
#Not working: setting geometry and reprojecting for second time.
gdf = gdf.set_geometry('geom2')     #still in crs_lonlat...
gdf.crs                             #...but this still says crs_new...
Out: 'epsg:3395'

gdf = gdf.to_crs(crs_new)           #...so this doesn't do anything! (geom2 unchanged)
gdf
Out:
                             geom1                      geom2
0  POINT (1001875.417 6948849.385)  POINT (8.00000 63.00000)
1  POINT (1001875.417 7135562.568)  POINT (8.00000 64.00000)

Ok, so, apparently, the .crs attribute of the gdf is not reset to its original value when changing the column that serves as the geometry - it seems, the crs is not stored for the individual columns. If that is the case, the only way I see to use reprojection with this dataframe, is to backtrack: start --> select column as geometry --> reproject gdf to crs_new --> use/visualize/... --> reproject gdf back to crs_lonlat --> goto start. This is not usable if I want to visualise both columns in one figure.
Second try
My second attempt was, to store the crs with each column separately, by changing the corresponding lines in the script above to:
gdf = gpd.GeoDataFrame()
gdf['geom1'] = gpd.GeoSeries([Point(9,53), Point(9,54)], crs=crs_lonlat)
gdf['geom2'] = gpd.GeoSeries([Point(8,63), Point(8,64)], crs=crs_lonlat)

However, it soon became clear that, though initialised as a GeoSeries, these columns are normal pandas Series, and don't have a .crs attribute the same way GeoSeries do:
gdf['geom1'].crs
AttributeError: 'Series' object has no attribute 'crs'

s = gpd.GeoSeries([Point(9,53), Point(9,54)], crs=crs_lonlat)
s.crs
Out: 'epsg:4326'

Is there something I'm missing here?
Is the only solution, to decide on the 'final' crs beforehand - and do all the reprojecting before adding the columns? Like so...
gdf = gpd.GeoDataFrame(crs=crs_new)
gdf['geom1'] = gpd.GeoSeries([Point(9,53), Point(9,54)], crs=crs_lonlat).to_crs(crs_new)
gdf['geom2'] = gpd.GeoSeries([Point(8,63), Point(8,64)], crs=crs_lonlat).to_crs(crs_new)
#no more reprojecting done/necessary/possible! :/

...and then, when another crs is needed, rebuild the entire gdf from scratch? That can't be the way this was intended to be used.
(Full disclosure: I posted this question on stack overflow, but it wasn't getting much love/attention over there.)
EDIT: Answer - or workaround
geopandas does not accommodate this use case at the moment. My workaround is to not use a GeoDataFrame, but rather combine a normal pandas DataFrame, for the non-shapely data, with several geopandas GeoSeries, for the shapely geometry data. The GeoSeries each have their own crs and can be effortlessly reprojected whenever necessary.


texts: ogr2ogr export to KML from SQLite, when using SELECT for fields Linestrings are lost

texts:
I'm using:

ogr2ogr -sql "SELECT * FROM 'trails_utm_utf8' ORDER BY name ASC" -f KML trails.kml trails_utm_utf8.sqlite -s_srs EPSG:900913 -t_srs EPSG:4326

to reorder, re-project and export to KML from a vectorlayer with linestrings kept in a SQlite.db. Now, I'd like to ommit the creation of the ExtendedData by the KML driver and tried so by putting the columns 'name' and 'description' from my SQLite file to SELECT - but strangely this will end up in a KML without any linestrings from the input file.

ogr2ogr -sql "SELECT name, description FROM 'trails_utm_utf8' ORDER BY name ASC" -f KML trails.kml trails_utm_utf8.sqlite -s_srs EPSG:900913 -t_srs EPSG:4326

Any ideas how to export only the fields/columns 'name' and 'description' and avoid that all other fields from the input file are put to the ExtendeData?


texts: Lable line points with geographical coordinates in QGIS

texts: I have a .shp file set of Congressional districts as lines.  Since the lines are composed of geographic coordinates, is it possible to label a point on a line and have the label automatically show the coordinates?


texts: How to get an elevation profile for a gps track?

texts:
I'd like to get a reasonably accurate elevation profile for a track recorded with a GPS (which often has very unreliable altitude data and occasionally none at all, depending on the model.)

Does anyone have any hints on the easiest way to do this. The two techniques I am considering so far are:

Using the Google Elevation API

This API is relatively easy to use, but still requires a few steps that aren't trivial due to it's usage restrictions: max 512 samples returned per request, and the number of points along the path is limited (by URL length) as well.

I expect a gpsbabel simplify filter can be concocted to reduce the track to a suitable number of points (no point in them being closer than 100m or so together due to the resolution of the altitude data), but then the problem remains of how to map this simplified track back onto the original path, since the lengths will differ.

Or, if this isn't suitable for automation the best approach may be to let the user select the transect points on a map manually.
Downloading the Shuttle Radar Topography Mission (SRTM) data and doing the query locally.

This is something I have no experience with, so any suggestions on how feasible this is are welcome.  How big is the data set? What GIS software is required, and can it be scripted in a suitable manner? I'd prefer not to have to write an sampling and interpolation algorithm, that sounds like a pain. What is the likely performance of such an approach? (I need it to be pretty quick and run on a memory-limited VPS webserver...)

Some further details to flesh out @MerseyViking's answer re downloading the data from http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp:

There are 72 x 24 tiles, each about a 20mb zip file that decompresses to a 72.1mb 16bit TIF file (6001x6001 pixels).

That's ~120 gb, which is more than I can store. Leaving it compressed and ignoring the oceans will reduce it to maybe 10gb, which is still a bit too large. Loading the data on-demand would dramatically reduce the needed storage space, but the source site is slow (I was only getting 10kb/s) making that pretty impractical.


texts: What are the PostGIS Raster Constraints used for?

texts:
My question is quite basic, but I have not found a simple explanation of why raster constraints are useful.

I have worked with rasters in PostGIS for a while but only recently I realized that I do not understand what the raster constraints are for. I understand that the constraints may provide some general information about the rasters in a specific table, but I have never accessed the constraints directly and am wondering how they are used internally in PostGIS.

My motivation for asking this question is that I have a single raster column with rasters of different sizes, extents and SRIDs and so the raster constraints are not very well bounded. If the raster constraints are used for indexing this might have implications on performance of the raster table in terms of speed when reading from the table.


texts: Works in Console but not in Script Editor

texts:
On QGIS, I am using a simple script to zoom the selected features. I write them in the Text Editor, but it didn't work until I used the console... and it does work. What kind of problem is this? Is there a token or anything I should introduce in the text editor?

municipios = QgsVectorLayer("\\","name","ogr")
QgsMapLayerRegistry.instance().addMapLayer(municipios)
l = iface.mapCanvas().layers()
expr = QgsExpression( "" )
it = l.getFeatures( QgsFeatureRequest( expr ))
itera = [i.id() for i in it]
l.setSelectedFeatures( itera )
iface.mapCanvas().zoomToSelected()

This is the error message (but not having it in my console, only with Text Editor of QGIS)

  Traceback (most recent call last):   File "", line 1, in
     File "c:/users/miguel~1/appdata/local/temp/tmpmlvs3w.py",
  line 9, in
      it = l.getFeatures( QgsFeatureRequest( expr )) AttributeError: 'list' object has no attribute 'getFeatures'


texts:
I want to load CSV file to QGIS. It works if I use the strings below:
uri ="file:///C:/Users/User/Desktop/Alex working/Tracks//Track.csv?delimiter={}&crs=epsg:25831&xField={}&yField={}".format(",", "Easting", "Northing")
vlayer = QgsVectorLayer(uri, "layer name you like", "delimitedtext")
if not vlayer.isValid():
    print("Layer failed to load!")

But if I pass Filepath and name of the file to variable, the file fails to load.
import os
from qgis.core import (QgsVectorLayer)
Path="C:/Users/User/Desktop/Alex working/Tracks"
name="Track.csv"
uri ="file:///Path/name?delimiter={}&crs=epsg:25831&xField={}&yField={}".format(",", "Easting", "Northing")
vlayer = QgsVectorLayer(uri, "layer name you like", "delimitedtext")
if not vlayer.isValid():
    print("Layer failed to load!")

What am I missing?


texts: Move Nominatim data from one system to another system

texts:
We have a Ubantu 14 server which is running the Nominatim service with postgres database.

Now we are tring to upgrade the Ubantu to latest, not able to upgrade because the postgres is installed separate, not used the "apt-get" which is not allowing to upgrade.

So we are planning to create new VM with Ubantu 16 and install postgres from apt-get.

Question is can i move the Nominatim database to new server? and How?


texts:
I was given a list of addresses without any zip codes, and I'm trying to write a python snippet that will loop through a set of possible zip codes and make a set of geocoded shapefiles that I will later merge.

My address table is essentially laid out like this: ("Project_" fields are the possible zip codes)

Address          City        Project_1   Project_2  ...  Project_42
123 Main Street  Birmingham  78954       78953      ...  78961

My python snippet looks like this:

import arcpy
from arcpy import env
env.workspace = "N:\Reference\GIS\State\Project\Address Locator"
geocoder = env.workspace + "\State_Addresses_Create_Addr"
for y in range(1,43):
    zippo = "Projects_" + str(y)
    outfile = '"N:\Project_Geocode_Results_'+ str(y) + '"'
    fncgeo = '"Street Apartment VISIBLE NONE; ZIP ' + zippo + ' VISIBLE"'
    arcpy.GeocodeAddresses_geocoding("Geocode_Table", geocoder, fncgeo, outfile)

When I try to run the snippet outside of a loop (i.e., based on a single field), the geocode works fine. When I run it as written above, I get:

  ERROR 000010: Geocode addresses failed. Failed to execute (Geocode Addresses).

Does anyone have any ideas as to why the geocode won't work in a loop?


texts: Separate a string value in two integers in QGIS

texts:
My problem is the following.
I have a huge dataset of forest stands.
All the information are separated but in many cases the stands have different age and this is marked in one attribute
for example

age = 14/2

This marks that one part of the stand is 140 years old and the rest 20 years.

My problem is that I need to have these two values separate, optimally in two colums:

age 1 = 14
age 2 = 2

Also, not always is the first number consisted of two digits, and vice versa.

Any ideas how I can do it?
I was thinking of using FieldPyCalculator, but I'm not any good with Python.


texts: "shutil.rmtree" doesn't get rid of my intermediate folders

texts:
I have a script that works fine, but now I'd like to remove all the intermediate folders I created therein in order to save memory and to be clean.
I've been trying using the "shutil.rmtree ()" function, as I've always done before, but it seems not working this time. Why?
Here is my script (the function is at the very end):

import arcpy, os, sys, shutil

try:
    newpath_01= r'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_01'
    if not os.path.exists (newpath_01):
    os.makedirs (newpath_01)

arcpy.env.workspace = 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\MODIS_Clip\\Sample_test'
    arcpy.env.overwriteOutput= True
    from arcpy.sa import *

    listMODIS = arcpy.ListRasters ()
    arcpy.CheckOutExtension ('Spatial')

    for MODIS in listMODIS:
        outSnowmask = ExtractByAttributes(MODIS, "VALUE = 200")
        outSnowmask.save (newpath_01+ '//'+ 'snow_'+ MODIS)

    newpath_02= r'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_02'
    if not os.path.exists (newpath_02):
        os.makedirs (newpath_02)

    arcpy.ClearEnvironment ('workspace')
    arcpy.env.workspace = newpath_01

    listMask = arcpy.ListRasters ()
    aspect_reclass= 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\DEM\\per_MODIS\\aspect\\asp_rcls_cknp'

    for mask in listMask:
        arcpy.CheckOutExtension ('Spatial')
        outExtractByMask = ExtractByMask (aspect_reclass, mask)
        outExtractByMask.save (newpath_02+ '\\'+ mask)

    newpath_03= r'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\temporary_03'
    if not os.path.exists (newpath_03):
        os.makedirs (newpath_03)

    arcpy.ClearEnvironment ('workspace')
    arcpy.env.workspace = newpath_02
    arcpy.env.overwriteOutput= True

    listExtraction= arcpy.ListRasters ()
    elev_range= 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\DEM\\per_MODIS\\recls_CKNPdem.lyr'

    for extraction in listExtraction:
        arcpy.CheckOutExtension ('Spatial')
        name_MODIS= str (extraction [18:25])
        ZonalHistogram (extraction, 'ASPECT', elev_range, newpath_03+ '\\'+ name_MODIS+ '.dbf')

    arcpy.ClearEnvironment ('workspace')
    arcpy.env.workspace = newpath_03

    listTable = arcpy.ListTables ()
    print listTable
    out_dir= 'C:\\Users\\unimi\\Documents\\Umberto\\Universita\\PhD\\Karakoram\\MODIS_Pakistan\\Aspect\\test_CKNP_aspect\\'

    for table in listTable:
       arcpy.AddField_management (table, 'FILENAME', 'TEXT')
       arcpy.CalculateField_management (table, 'FILENAME', '"'+ table+ '"')

    arcpy.Merge_management (listTable, out_dir+ 'Snow_cknp_sample.dbf')

    shutil.rmtree (newpath_01, newpath_02, newpath_03)

except:
    print arcpy.GetMessages ()

Any help will be very appreciated.


texts: Preventing ArcMap layout view re-drawing every time label moved etc?

texts:
In the Layout View whenever I move graphics, or labels within it, ArcMap  10.1 redraws it.

How do I turn that off to prevent it from re-drawing repeatedly?


texts: Getting NDVI values using Google Earth Engine Python API

texts:
I am new in GIS. I want a function in Python API that will give me NDVI values of a given region within specified dates. I have tried the following code but print returns None.
dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_8DAY_NDVI').filterBounds(ee.Geometry.Point(-122.1, 37.2)).filterDate('2019-01-01', '2019-12-31')
colorized = dataset.get('NDVI').getInfo()

print(colorized)


texts: I'm using QGIS 3.8 (the latest version) and I want to reproject a shapefile from WGS84 to North America Conformal Conic. This shapefile contains locations of the stations (points). However, I spent many times to do that, but it does not work because the reprojected shapefile showed me the stations in a single lines. In ArcMap, it was easier. Maybe I missed something.


texts: How to determine the average NDVI value of a selected area?

texts:
I'm new to ArcMap 10.2 and trying to get NDVI data of city centers.

So my problem is, after I'm calculating the NDVI and getting the result raster, I would like to get a mean (average) value of specifically and only the area I select (which is the city center on a 1:50000 scale).

Things I tried and didn't make me happy;

1) Identify tool. It constantly gives me the information about only ONE pixel even if I draw an area on raster image with shift + click.
2) Clipping raster from View > Data Frame Properties > Data Frame > Clipping tool. Even after I crop the satellite image to 1:50000 area I'm interested in, I can still get information from blank places OUTSIDE of my area, and it also affects the mean value of AOI NDVI.

So basically, I want to get an accurate mean NDVI value ONLY including the pixels of the selected area of interest.

I'm using LANDSAT 7 and 8 satellite images which are way larger than my AOI.

Sorry for bad English.


texts: XML entities not working with node-mapnik

texts:
I'm not getting XML entities to work with node-mapnik. In the changelog for 1.4.4 I read:

"Mapnik now build with rapidxml/ptree xml parser rather than libxml2 (meaning no xincludes or entities are supported)".

Are there other ways of using variables or include other XML stylesheets in node-mapnik? Is it possible to configure node-mapnik to support XML entities and includes?


texts: How to add Simple 3d Bar Graph on top of 2d Feature class in ArcMap

texts:
I am trying to add a simple bar graph on top of a map in ArcMap. The idea is very simple, to calculate x,y (suppose building area vs building name) and create a bar graph of it. This can easily be done using create graph tool but I want to show a 3D graph(vertical Bars) on top of flat 2D shapefile/feature class. I tried arcScene but this feature is not activated in that. I can create bar graph in 'Symbology' tab but in Arcmap but it won't let me have a 90 degree angle between the graph and feature class.
Is it possible in ESRI Arcmap with or without 3D analyst?

Example of a map I want


texts: Tree Cover Change analysis on QGIS, "Stripes" or "Banding" on the image

texts:
I am currently looking at Tree Cover Change around a set of 288 villages in the Equatorial Forest of Africa using QGIS.

In order to do so I am using the GFCC30TC: Global Forest Cover Change (GFCC) Tree Cover Multi-Year Global 30 m V003 dataset. This dataset contains data for the years 2000, 2005, 2010 and 2015.
Please find information on the data here:

https://lpdaac.usgs.gov/dataset_discovery/measures/measures_products_table/gfcc30tc_v003

I am basically using a Difference Image Approach, subtracting images band pairs for each year (2000-2005), (2005-2010) and so on to highlight changes.

However, the results I find are biased because the images for 2005, 2010 and 2015 contain "stripes" or "banding" leading to biased calculations. Here is an example in the picture below (raster of change between 2000 and 2005).

The problem for me is that pixels in the stripes take the value of 0 and thus bias my analysis.

Does anyone have any solution to overcome this issue?

The best for me would be either to vanish the stripes or to find a solution to set the stipes value to missing and to account for this in the calculation.

As I am not a GIS expert (by far), maybe some of you know any other satellite imagery database that do not have such stripe issue.


texts: How to determine if lat/long is within SABS school boundary?

texts:
I am an absolute beginner with GIS data.   I am attempting to write a tool to determine if a lat/long point resides within a polygon within a shapefile provided by SABS (https://nces.ed.gov/programs/sabs/).  Basically, I am trying to write something that will tell me the primary school for a given address when internet is not available (no Google Maps, GreatSchools, etc.).

Anyway, I am using GeoTools with modified PointInPolygon code from https://github.com/ianturton/geotools-cookbook/blob/master/modules/filters/src/main/java/org/ianturton/cookbook/filters/PointInPolygon.java to search for lat/long points in polygons.

With that, I can see that there are lat/long points within the GeoTools suggested datafile (http://www.naturalearthdata.com/downloads/50m-cultural-vectors/).  However, the SABS data does not report any points in polygons.

When I look at the polygon data for the two files, the points for the suggested data file are in the lat/long ranges I expect (+/-90, +/-180).  The point data in the SABS file has values many magnitudes larger.

My guess is that the SABS data is in a projection that doesn't play well with querying if a point is inside the polygon.  I don't know how to transform the data so that I can perform lat/long queries.

Can someone explain what I need to do to be able to perform a lat/long search with the SABS data?

SABS prj: PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Albers"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-96.0],PARAMETER["Standard_Parallel_1",29.5],PARAMETER["Standard_Parallel_2",45.5],PARAMETER["Latitude_Of_Origin",23.0],UNIT["Meter",1.0]]

Natural Earth prj:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Code added to PointInPolygon to test data files:

for (double lat = -90.; lat < 90.; lat += 0.01) {
        for (double lon = -180; lon < 180; lon += 0.01) {
            final Point p = fac.createPoint(new Coordinate(lat, lon));
            final boolean flag = tester.isInShape(p);
            if (flag) {
                System.out.println(p + " is inside");
            }
        }
    }


texts:
I can't get my head around plotting a circle on a big sphere (Planet Earth).
I have a column that has 180 rows - 1 for each degree.
The next row is each degree converted to radians.
Easy so far.

A1 is the centre point of the radius latitude

B1 is the centre point of the radius longitude

C1 is the radius required in miles

Column D is each degree so there are 180 rows.

Column E is each degree converted to radians

Column F has (first three rows):

A1+((C1/69)*COS(E1))    latitude for degree 1.
A1+((C1/69)*COS(E2))    latitude for degree 2.
A1+((C1/69)*COS(E3))    latitude for degree 3.

Column G is (first three rows):

B1+(C1/(COS(RADIANS(A1))*69)*SIN(E1))   longitude for degree 1.
B1+(C1/(COS(RADIANS(A1))*69)*SIN(E2))   longitude for degree 2.
B1+(C1/(COS(RADIANS(A1))*69)*SIN(E3))   longitude for degree 3.

After the first 180 degrees, the formula changes as you multiply by -SIN & -COS.

A1+((C1/69)*(-COS(E1)))   latitude for degree 181.
A1+((C1/69)*(-COS(E2)))   latitude for degree 182.
A1+((C1/69)*(-COS(E3)))   latitude for degree 183.

and

B1+(C1/(COS(RADIANS(A1))*69)*(-SIN(E1)))   longitude for degree 181.
B1+(C1/(COS(RADIANS(A1))*69)*(-SIN(E2)))   longitude for degree 182.
B1+(C1/(COS(RADIANS(A1))*69)*(-SIN(E3)))   longitude for degree 183.

I get a wonderful circle in the graph which overlays a map, can be moved and resized by changing A1, B1 & C1 (A1 & B1 use a filterxml-webservice, and it's actually just a postcode I change), and looks very good.

My problem is, have I actually got the math right?

Is it correct to divide the radius required in C1 (e.g. 10 miles) by the distance per degree at the equator (app 69 miles)?

I personally live around 53deg latitude, so for me 1 deg longitude is about 41.5 miles but the spreadsheet is used worldwide.

I've removed the $ signage for clarity (used in spreadsheets to drag formula down rows & across columns).

I've been asked to give a screenshot.  Don't know what good it'll do, as I'm just asking if my sums are right...

Two circles are shown with the cluster of members in each (the names are fake - just in case you were wondering)
I believe the circles (especially the lower one) are too big  The radius is set to 35 miles, but it's looking more like 45. If the circles were correct, there would be 15 more members on the sheet

I came to this site as it was here that I found some calculations to help me at: Calculating longitude length in miles?

The very simplified sheet is here


texts:
I'm trying to calculate EVI in Arcmap using raster calculator. I put in this calculation

2.5 * ((Float("Band4") - "Band3")) / (Float("Band4") + 6 * "Band3" - 7.5 *"Band1" + 1))

and the output is syntax error. Do you know what I'm doing wrong?


texts:
I'm a student studying GIS, and I need to create geocoding for my work.

How can I create the address locator style for fields with my requirements and my specified shapefile.

My attributes are:

Id
Survey_No
Owner_Name
Cadestre_No
Soil_Type
Cult_Type
Water_Sour
Rights
Bebts


texts: Difference between feature to raster and polygon to raster

texts:
(I'm using Arc 10.1)

I'm trying to rasterize a data set. I have a shapefile with geographic coordinates (GCS_GDA_1994) which I want to convert to a raster with 100m grid cells. So I need to project into a projected SRS - I used GDA_1994_South_Australia_Lambert but I don't think it matter which one for my question. So I did that, and afterwards I tried to use 'Polygon to raster' to convert to a raster. However, in the input field 'Cellsize' I get by default a very small value (0.0002 or so) - reading the help file, this makes me thing it is because it is still using degrees as the unit there. And indeed when I put in 100 there (for 100 meters), I get one big cell; it is using the wrong unit. When I check 'linear unit' under the layer properties, it does say 'Meter' however.

So then I tried Feature to Raster and lo and behold, I get a default suggestion in the correct order of magnitude and putting in 100 manually gives me a 100m raster.

Now, Polygon to raster has more options that control how the value of the cell is calculated, and I need that. So my question is: what is the difference between these two tools that makes one use a different unit than the other, and how can I make 'Polygon to raster' use meters?


texts: In ArcMap, when I use PointsToLine in ArcToolbox, error appears

texts:
In ArcMap, I open the Python window and print some codes as followed:

import arcpy
arcpy.env.workspace = "C:/data"
arcpy.PointsToLine_management("calibration_points.shp",
                              "C:/output/output.gdb/out_lines",
                              "ROUTE1", "MEASURE")

  Runtime error   Traceback (most recent call last):

  File "", line 1, in

  File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy__init__.py", line 21, in

      from arcpy.geoprocessing import gp

  File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing__init__.py", line
  14, in
      from _base import *

  File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing_base.py", line
  598, in
      env = GPEnvironments(gp)

  File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing_base.py", line
  595, in GPEnvironments
      return GPEnvironment(geoprocessor)

  File "c:\program files
  (x86)\arcgis\desktop10.3\arcpy\arcpy\geoprocessing_base.py", line
  548, in init
      import weakref

  File "C:\Users\yanglan\Anaconda2\Lib\weakref.py", line 14, in

      from _weakref import ( ImportError: cannot import name _remove_dead_weakref

i want to know whether this error has a correlation with the anaconda2 installed in my computer, so curious about how to fix this problem


texts: Calculating distance between borders when one polygon inside another polygon?

texts: I am brand new to ArcGIS and having some trouble. I would like to calculate the distances between the border of one polygon to another polygon. However, these polygons overlap.  Thus when I use the Near function, it results in zero because the polygons overlap at some points or one is within a larger one.  I cannot use the Point Distance because I have shapefiles and so I have nothing to choose from in the drop down menu.


texts: Add features to main SDE Table with geotools

texts:
Sorry for my english, I hope you can understand ;)
I'am using geotools 8.4 version. I try to copy all features from SDE Table(TEST.DBO.MM_SOURCE) to another SDE Table(TEST.DBO.MM_TEST). All features are copied correctly, but they are not stored in main table (TEST.DBO.MM_TEST). When I open ArcMap, edit some feature in TEST.DBO.MM_TEST and save, all features are back in main table. Both SDE tables are registered as versioned and moving edits to base. What I must do to store all features in main table?

Here is my java code:

String server = "localhost";
String port = "5151";
String database = "Test";
String user = "sa";
String password = "sa";
String version = "dbo.Default";

connector = new SDEConnector();
connector.setConnectionParameters(server, port, database, user, password, version);
connector.connect();

SimpleFeatureStore dataStore =
        (SimpleFeatureStore) connector.getDataStore().getFeatureSource("TEST.DBO.MM_TEST");

SimpleFeatureCollection fc1 = connector.getFeatureCollection("TEST.DBO.MM_SOURCE")
SimpleFeatureCollection fc2 = dataStore.getFeatures();

SimpleFeatureIterator iterator = fc1.features();
while (iterator.hasNext()) {
    SimpleFeature feature = iterator.next();
    Geometry geometry = (Geometry) feature.getDefaultGeometry();

    SimpleFeatureBuilder builder = new SimpleFeatureBuilder(fc2.getSchema());
    builder.set("Shape", geometry);
    builder.set("P1", feature.getAttribute("P1"));
    builder.set("P2", feature.getAttribute("P2"));
    builder.set("Source_OID", 829023);
    builder.set("Marker", 10);

    SimpleFeature sf = builder.buildFeature(null);
    // SimpleFeature sf = builder.buildFeature(builder.createDefaultFeatureId());

    List<SimpleFeature> list = new ArrayList<SimpleFeature>();
    list.add(sf);

    dataStore.addFeatures(DataUtilities.collection(list));
}
iterator.close();
time = System.currentTimeMillis() - time;
System.out.println(time);


texts: Extracting value from three different raster for same pixel/position as CSV or tabulated dataset using ArcMap

texts:
I have computed Hydrologic Sensitivity Index (HSI) as a raster for a region using slope, Hydraulic conductivity, Depth to the Restrictive Layer and Flow accumulation. All the files are in rasters of same pixel size. The rasters are float type rasters.
Now, I would like to plot Hydraulic conductivity vs HSI, Depth vs HSI for every pixel etc.. In other words, for each pixel, I need the computed HSI value, Depth, Hydraulic conductivity, and Flow accumulation in a column or dataframe.
How can I accomplish this in ArcMap?
What I tried:

First I tried to use the Combine tool. But Combine tool doesn't seem to work with float raster even if I convert it into integers there will be too many unique values to make this combination.

What if two different pixels had same values?
In that case, I do not think this tool will show those two pixel separately (right?).

Convert each Raster to Int Raster than to Polygon than use Union tool. I couldn't accomplish this because the HSI raster file is very detailed and takes too long to convert to polygon and seems to crash in the middle.

What I think could work:

I know there is a way to convert Raster data into a Numpytable. This way I will have four separate Numpy table for each of the four data but how will I map the value in such a way that i ensure that they are for the same pixel?

I know I can use arcpy.Describe.polygon.OIDFieldName to create id in Polygon but is there similar tool for Raster. But i am also not sure when I am using this tool for four different raster.
Will the pixel in the same location get the same oid values?


texts: Leaflet and census block data: not displaying

texts:
I am new to leaflet and am trying to load a geojson file showing census blocks and associated information for a city in CA. It is a small city and not too many census blocks. I cannot get the popup feature to display any data. Values are written to the console but not showing up in the map. I get a warning here:
Invalid CSS property value: transform 0.25s cubic-bezier(0,0,0.25,1)

My code is:

<!DOCTYPE html>
<meta charset="utf-8">

<link rel="stylesheet" href="leaflet-0.6.4/leaflet.css" />


texts:
In layout view I have many pages with different parts of the main map. For each of these I have a legend.
Is it possible to designate all of them the same kind, and when I update one I update all?


texts:
I am new at open source GIS platform. I am using Ubuntu 12.10. I have tried to install the followings but without 2-3 tools I could not install anyone. Can any body help me how can I install the following software from scratch. I mean dependencies on each other and installation priorities.

PostgreSQL
PostGIS
GDAL
GEOS
PROJ4
PGRouting (cgal & gaul)
Tilemill
Mapnik
Mod_tiles
QGIS


texts: KML Regions feature not working in Google Earth 8.0.2.2334 installed on Android

texts:
When I load a kml in Google Earth 8.0.2.2334 installed on my Android phone, the KML Regions feature is not working (for me).
I have read and reread the Google Dev pages. I have googled for a solution for two days.

I am really hoping somebody will simply alter one of the examples below and then I will load it to my phone and see that it is possible to set limitations (zoom dependent) to when my data (kml) is visible in Google Earth on my phone. or load an example of their own.

I want to keep this simple so I will use existing samples to show how the feature  works as intended on my and most likely your PC version of Google Earth but then not at all on my phone's Google Earth app.

Sample Examples can be found at Google's KML Samples webpage, just scroll down that page to find the "Region KML" examples

Here are two of them:

http://kml-samples.googlecode.com/svn/trunk/kml/Region/polygon-simple.kml

http://kml-samples.googlecode.com/svn/trunk/kml/Region/polygon-swap-pop.kml

As you might have read on Google Developers website:

"Regions are a powerful KML feature that allows you to add very large datasets to Google Earth without sacrificing performance. Data is loaded and drawn only when it falls within the user's view and occupies a certain portion of the screen. Using Regions, you can supply separate levels of detail for the data, so that fine details are loaded only when the data fills a portion of the screen that is large enough for the details to be visible."

I want to know how to make them work on my phone please.
In case you ask:
Android ver. 5.0.1
HTC Sence ver 6.0
Google Earth 8.0.2.2334


texts:
Is there any way to report all the steps I did within a ArcMap 10 project. this could help me to:

Make a Metadata description
Trackback possible errors afterwards and/or control the analysis
Easier repeat the analysis some other day

When working with scripts you don't need that beacause they clearly describe what has been done but when clicking around it would be fine to be able to see afterwards what was done. So instead of writing it all down to a textdocument I'm searching for an easy and convinient way to save this automatically.


texts: How can I implement Persistance in a Layer Extension?

texts:
Last time I checked, arcgis doesn't honor IPersistVariant in layer extensions.  Rather advanced code has been provided here showing how to implement IPersistStream.

Is there something simpler?


texts:
Is there any public or private sources of LIDAR data for the city of Chicago? Searching online has not lead me to any sources, and I am unsure how the commercial world of LIDAR works.

The LIDAR point spacing would need to be at an interval that I can use to estimate the building height from the footprints.


texts: Leaflet map isn't loading in the correct location

texts:
I'm currently working with Leaflet and a dockable web framework called Golden Layout. I'm trying to load tiles into the map, but they are not loading properly and I'm not sure why.

Currently, I have created a function to handle the creation of the map. The map container itself is created correctly, however when I try to add the tile layer, that's where things go wrong.

There are a couple problems with adding the tile layer. The first is that the tiles don't load except for the one in the top left hand corner. The second problem is what tile is being loaded. The tile that gets loaded in the top left hand corner is the tile I have set to be the center of my map.

If the page is adjusted in anyway (ex: console opened), all the tiles will load and the center will be correct.

I've created a Codepen to show the problem.


texts: Save QGIS point layer as KML and retain attribute table

texts:
QGIS 2.18.5

I've done this in the past but I can't remember how I did it! I'm saving a point layer as a kml and trying to bring over the attribute data across with the kml. I've got some screenshots of the result from the past:

I've read a post about exporting 2 attributes but obviously I managed to export a selection of my attributes. I'm not an experienced coder and there weren't any intricate expressions involved. I remember filling in the 'Description Field' and 'NameField' shown below but from recollection, that only changed what the pins in Google Earth are titled with - shown as parcel numbers in my case.


texts: Get current zoom level from QGIS map canvas

texts:
Is there any particular way of getting the current zoom level from the map canvas through python code in QGIS environment?
I want to convert coordinates of the cursor and the current zoom level into Google Map Tiles coordinates.


texts: Adding GeoRSS feed as point layer using OpenLayers?

texts:
There may be a simple answer, but I have been looking on the internet for answers and trying different things.

I want to display this RSS feed as a point layer in OpenLayers 2.11.

Can anyone please give me some code that may help?


texts: Batch export multiple datasets in a group within a (mesh/raster) layer as images?

texts:
I have a mesh layer that contains a group with multiple datasets (dates). I'd like to export the visualization of each date together with another active raster layer (showing the basemap).
I could do this manually by selecting the dataset and then click on "Export Map to Image". This would take me ages though.
Any suggestions?


texts: Boundless Suite - Composer Security Permissions

texts:
I’m having trouble getting users-roles-data access rules correctly configured in the Boundless Suite 4.9.1.  The issue appears to be with Composer, not Geoserver.

BLUF:  I want 2 separate editors for 2 separate workspaces – I don’t want either to be able to upload data, view, make maps, edit or delete data in the other workspace. I can’t prevent a user from viewing or uploading data to a workspace.

Simple use case, in Composer I want:

One user to be able to edit a private workspace, but not view or edit data in the public workspace.
One user who is able to edit a public workspace, but not view or edit data in the private workspace.

I am testing behavior by seeing if users can upload data and view data in the workspaces based on their permissions. My goals is that the public user would not be able to upload data to the private workspace and vice/versa. I started with the following data access rules:

There is a 1:1 correspondance between the privateEditor user and the c1_role_privateEdit and between the publicEditor user and the c1_role_publicEdit  In this case both users can upload, view and delete data in both workspaces. I assume the access rules are where the problem is but I have tried every configuration I can think of:

1.      I have locked the write access down for the admin user but no effect.

(..R=*; ..w=admin,groupAdmin);

2. Adding read permission to any user to any workspace results in no longer being able to upload data to that workspace, even if they also have write permission. For example I add read permission to the private user in the private workspace (C1_private..w, C1_private..r) – which the private user is not able to upload data to the private workspace but the public user can still upload data to the private workspace. When I restrict the Admin user to read access (..r=admin,groupAdmin) due to inheritance the other users are unable to upload data to any workspace.

I have added additional admin privileges to the private user (C1_private.*.a) for the private workspace, and the public user for the public workspace with no effect.
I have tried making he role_publicEdit the parent to role_privateEdit with no effect.
I’ve tried locking services under the unlikely hypothesis that uploading and deleting may call WFS or WPS but service level permissions are granted on all workspaces and I want to restrict specific users to specific workspaces.
I looked at the application calls to see if I needed to lock down permissions in the rest.properties XML but as far as I can tell uploading and downloading are not running through the rest api. I have tried setting up every combination of r, w, and a for users to those workspaces as the admin user in Geoserver but no matter what I do, within Composer each editor can see, upload to, and delete from workspaces that are not assigned to them.

Is there another setting that I need to be looking at in order to lock down the security in Composer? Are uploads and delete considered write functions, or admin functions?


texts: QGIS merge vector layers and add a group field

texts:
Is there any way to merge several vector layers into a single layer while adding an attribute column that groups the points by the layers they belonged to prior?

So I have multiple delimited txt files with coordinates loaded into my project. I have figured out a couple of ways to merge them into a single shapefile but a group # or name would make it a lot quicker for me to use Points2One to create line data from my many sets of points.


texts:
Is it possible to perform a vector machine classification of raster image using the ArcGIS software package. This supervised classification method is available in the ENVI software but not for ArcGIS.

I have integrated the ENVI toolbox into my ArcCatalog Toolbox set. While the ENVI toolbox offers me the tools to perform supervised classification using four different methods, it does not include the Vector Machine Classification Method.

Since, I've access to an ArcGIS license but not an ENVI license, I am looking to access this functionality in ArcGIS. Any suggestions on how this can be achieved?


texts: Projected CRS interfering with Spatial Analyses?

texts:
I am working with two polygon shapefiles, one of which is a buffer around the other, and a point shapefile in a project set to EPSG:54016. All layers are correctly reprojected and in the same CRS, OTF is turned off.
I am now trying to identify/select all points that lay within the buffer zone, but none of the tools work for me. Points in Polygon tool won't even start, Spatial Query and Select by location freeze the system.

I cannot find any information about it on the internet.

I wonder if anyone ever encountered this problem or if I am lacking knowledge about the procedures involved?


texts: Label moved with Move Label Tool disappear on QGIS v2.6.1

texts:
I have done a deep research on previous post about "move label tool", but I have not found the correct answer for my case.
I have the version 2.6.1 of QGIS.

When I select a particular point of my layer (because I need to move the label of this point), I can activate the "move label tool", but then when i try to move the label, the label disappear. I checked the attribute table and I have seen that the name of the label became as a coordinate. Do you know what I should do?
I have tried to use also the plug-in EasyCustomLabeling, but when I use it, the label  appears with a white buffer, that doesn't correspond to the buffer option of the label.

Any idea?


texts: Save wkt as "polygon_name".text from postgis table

texts:
I have a PostGIS-Table (PostGIS-layer) with more than 400 Polygons (ID, name, geom). I need the WKT of each feature in the layer saved in a "polygon_name".txt file with the name of the Polygon. How to solve that lovely by QGIS and pyQGIS or is there a direct way in postgis.

The result must be a folder with more than 400 *.txt files named by the names of the polygons.

After some help especialy by hannesW und jef I solved it in QGIS with pyqgis by the code:

    features = vlayer.getFeatures()
for feature in features:
    name = feature.attribute("name")
    filename = 'C:/temp/wkt/{x}.txt'.format(x=name)
    file = open(filename, "w")
    file.write(feature.geometry().asWkt())
file.close()


texts: How to calculate best path between 2 points from a DEM in Qgis?

texts:
I'm analyzing the roads and paths of a region by calculating how effective the communication is in the region.

First I calculated the shortest paths between the points using the actual road map.

Now I want to generate a new vector layer by calculating the best paths between the same points in the region from a digital elevation model, using parameters such as shortest distance, lower slope, etc..  The objective is to create artificial paths and compare them with the existing road net.


texts: Using standalone pyQGIS application QgsApplication.initQgis() fails?

texts:
I'm trying to execute the following code on a windows 7 system:

import qgis.core as qc
qgis_path = "C:\Program Files\QGIS Lyon\apps\qgis"
qc.QgsApplication.setPrefixPath(qgis_path, True)
qc.QgsApplication.initQgis()
print "Successfully Initialised"

The qc.QgsApplication.initQgis() command causes python to stop working.  Can anyone help with the steps I need to diagnose this problem?

Some background:

I've tried a number of different combinations of forward slashes, escaped backslashes etc.for the path and still get the same error;
I've been developing the code (part of a larger application) on Linux, without problems;
A long time ago, I tested the code on Windows and it worked OK.  At that time I had QGIS Dufour installed;
I suspect it may have something to do with having the correct environment variables set, but have been unable to find a definitive guide on which and how these variables need to be set.


texts: KeyError when adding layer to WebMap with add_layer in Python ArcGIS API

texts:
I'm a GIS student and in my school project I am trying to add a layer to a Web Map.  I am using Python ArcGIS API in Jupyter.
I found the layers I wish to add with content.search and have them in a list of tuples along with their titles.  Then I have this code to add the layers:

for item, title in featurelist:
    print(item)
    wm.add_layer(item, options = {"title":title})

It prints the item and then the Key Error:

<Item title:"Boat Ramps of Lake Travis by CR" type:Shapefile owner:cynthia.roush@g.austincc.edu>
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\gis\__init__.py in __getitem__(self, k)
   7004         try:
-> 7005             return dict.__getitem__(self, k)
   7006         except KeyError:

KeyError: 'layers'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
<ipython-input-87-120a1e45a0bb> in <module>
      1 for item, title in featurelist:
      2     print(item)
----> 3     wm.add_layer(item, options = {"title":title})

C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\mapping\_types.py in add_layer(self, layer, options)
    318             if not self._extent:
    319                 self._extent = layer.extent
--> 320             if hasattr(layer, 'layers'):
    321                 if layer.type == 'Feature Collection':
    322                     options['serviceItemId'] = layer.itemid

C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\gis\__init__.py in __getattribute__(self, name)
   6976     def __getattribute__ (self, name):
   6977         if name == 'layers':
-> 6978             if self['layers'] == None or self['layers'] == []:
   6979                 try:
   6980                     with _DisableLogger():

C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3\lib\site-packages\arcgis\gis\__init__.py in __getitem__(self, k)
   7007             if not self._hydrated and not k.startswith('_'):
   7008                 self._hydrate()
-> 7009             return dict.__getitem__(self, k)
   7010     #----------------------------------------------------------------------
   7011     @property

KeyError: 'layers'


texts: Copying exact field properties to new layer using ArcPy?

texts:
I fear I am missing something fundamental.

I can get a list of fields from one layer using ListFields, and look at the properties:

fields = arcpy.ListFields(feature_class)

# Iterate through the list of fields
for field in fields:
    # Print field properties
    print("Field:       {0}".format(field.name))
    print("Alias:       {0}".format(field.aliasName))
    print("Type:        {0}".format(field.type))
    print("Is Editable: {0}".format(field.editable))
    print("Required:    {0}".format(field.required))
    print("Scale:       {0}".format(field.scale))
    print("Precision:   {0}".format(field.precision))

Is there an easy way to add these fields with the exact same properties to a new layer using AddField_management?

I was hoping I could pass the field object directly and not have to decompose it into the individual properties. But the following does NOT work:

for field in fields:
    # add fields with the same properties
    arcpy.AddField_managment(feature_class2,field)

It does not seem very efficient to break out and pass every field property (field.name, field.type, etc). Surely there is a better way.

I get the error:

  ERROR 000623: Invalid value type for parameter field_name


texts:
I am reprojecting a SRTM 3 arc-second raster file from WGS84 geographic coordinate system, to projected coordinate system which is using Azimuthal Equidistant projection:

gdalwarp -s_srs EPSG:4326 -t_srs "+proj=merc +lat_ts=40.81266 +lon_0=14.414252" -r near -of GTiff C:/vesuvius2_wgs84.tif C:/vesuvius2_aeqd.tif

I attached the initial raster file in here (vesuvius_wgs84.tif).

The issue is that if I use the Nearest neighbor resampling method (-r near), some sort of "stripes" appear on terrain once I reproject the file to Azimuthal Equidistant projection. Here is the how the reprojected raster looks like:

And here is its 3d representation:

Interestingly similar "stripes" appear in both direction if I reproject the initial raster to UTM:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32633 -r near -of GTiff C:/vesuvius2_wgs84.tif C:/vesuvius2_epsg32633.tif

The "stripes" go away when I use Bilinear or Cubic resampling methods (-r bilinear, -r cubic) instead of the Nearest Neighbor one.

Why do these "stripes" appear when Nearest Neighbor resampling method is used?


texts: Is it possible to call a complicated query as a stored procedure within a geoprocessing service?

texts:
I need to implement a searching geoprocessing service.
A part of the query should be very complicated and use oracle text.
My service should support spatial search too.

1) I am thinking of performing the text query outside of arcgis by calling a sort of stored procedure from my toolbox. Is it possible ?
2) A solution could be  to construct a python tool that would assemble the sql quey using a range of parameters provided by the user. However I would like to avoid that solution since I dont want my query construction to take place inside the GP service.

Any ideas on what could be more efficient and easier or any alternatives?


texts: Missing files in Loader_Generate_Nation_Script for TIGER geocoder

texts:
I'm while running

CREATE SCHEMA
Archive:  tl_2015_us_state.zip
 extracting: /tmp/gis/gisdata/temp/tl_2015_us_state.cpg
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.dbf
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.prj
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.shp
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.shp.ea.iso.xml
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.shp.iso.xml
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.shp.xml
  inflating: /tmp/gis/gisdata/temp/tl_2015_us_state.shx
unzip:  cannot find or open */tl_*state.zip, */tl_*state.zip.zip or */tl_*state.zip.ZIP.

No zipfiles found.

This keeps happening for latter steps:

unzip:  cannot find or open */tl_*county.zip, */tl_*county.zip.zip or */tl_*county.zip.ZIP.

No zipfiles found.

Do I need to get those from a different source?

EDIT:

It seems to be a problem with the paths:

It downloads the files:

  ftp://ftp2.census.gov/geo/tiger/TIGER2015/ADDR/PLACE/tl_*_17_*
           => ‘ftp2.census.gov/geo/tiger/TIGER2015/ADDR/PLACE/.listing’
Resolving ftp2.census.gov (ftp2.census.gov)... 148.129.75.35, 2610:20:2010:a09:1000:0:9481:4b23

Can't go to the corresponding path:

 cd: can't cd to /tmp/gis/gisdata//ftp2.census.gov/geo/tiger/TIGER2015/ADDR/PLACE

The dir only exist until:

/tmp/gis/gisdata//ftp2.census.gov/geo/tiger/TIGER2015/ADDR/

It Doesn't find the downloaded files:

unzip:  cannot find or open tl_*_17*_place.zip, tl_*_17*_place.zip.zip or tl_*_17*_place.zip.ZIP.

Also, this other directory is missing:

No such directory ‘geo/tiger/TIGER2015/ADDR/PLACE’.

At which point should they be created?


texts:
In google map all the labels are rendering well and there are none missing. In my map some times only the point marker is showing and the label is not displaying.

I am asking about "Auto align and AUTO OFFSET" OPTION in geoserver like mapserver.

In this image all the labels are shown in different offset.

What do I have to change in my SLD?


texts: Seeking open source or free tool to open 3D Mesh GIS files?

texts:
I'm looking a way to open 3d Terrain Mesh files.

I have four different folders with same 3d content , I'd like to open one of them :

Cesium folder (with internal Tile folder, web.config, tileset.json)
OBJ folder (obj, obj.mtl and with internal Tile folder)
OSGB folder (with internal Tile folder)
XPL3 folder (with folders floor, Tile and LODTreeExport.xml)

I'm looking for open source or free tool to open 3d Mesh listed above.

I tried to open each of them in meshLab and BlenderGIS both don't have options to open these files are listed above, same above google earth,Global Mapper, openSCAD, QGIS, AutoCAD Map3D, Sketch Up and GRASS.

I am new to the field of GIS.


texts:
I am running a loop to go through tables of a database's feature classes to rewrite the URL links that are currently in the records. The first elif statement is to check whether there is an empty string or not and to proceed by doing nothing. I wrote this in because the first if statement that also accounts for NULL values does not check for empty strings. Instead, the script just uses the last elif and writes the URL with no filename at the end. All other conditions are fine.

parish = arcpy.GetParameterAsText(0)
arcpy.AddMessage('Parish Name =  --> "{}"'.format(parish))
field = 'Photo1'
path = r"http://website.gov/Photos/"

for fc in arcpy.ListFeatureClasses():
    arcpy.AddMessage('Processing Feature class -->"{}"'.format(fc))
    rows = arcpy.UpdateCursor(fc)

    for row in rows:
        photolink = row.getValue(field)
        arcpy.AddMessage('Current Photo1 value -->"{}"'.format(row.Photo1))

        if not photolink:
            pass

        elif photolink == "":
            pass

        elif photolink[:4] == "http":
            arcpy.AddMessage('Photolink If statement worked')
            photolinkMod = photolink.split(r'/')
            photo = str(photolinkMod[len(photolinkMod) - 1])

            arcpy.AddMessage('New Photo1 value -->"{}"'.format(photo))
            row.setValue(field, path + parish + r"/" + photo)
            rows.updateRow(row)

        elif photolink:
            row.setValue(field, path + parish + r"/" + photolink)
            rows.updateRow(row)


texts: Convert from Geometry to SDE Binary storage type in SDE database

texts:
We have been battling performance issues between SQL 2012 and ArcGIS 10.1. One of the most common solution is to convert to SDE Binary storage. We know that this storage type is not as dynamic as the Geometry type. But we just opted to try to convert atleast one table in the database (assuming that it is possible).

Problem is, we are not sure on the correct procedure to do this. Do we do it on the database side? or in the ArcGIS side?

Thanks.


texts: Getting number as parameter to use in SQL expression leads to ExecuteError: ERROR 000358?

texts:
In my script I am using the Select tool (arcpy.Select_analysis) and in the SQL expression, I want to use a number that I call in as a parameter.  I have trying to script it like:

subsize = arcpy.GetParameter(3)

where_clause = '\"AreaSqKm\" > {0}'.format(subsize)

arcpy.Select_analysis(waterbodies, lg_waterbodies, where_clause)

In this example I'm entering the parameter as the number 0.0005, as a double.

I've tried tweaking it various ways, but every time I try to run the tool I get:
ExecuteError: ERROR 000358" Invalid expression "AreaSqKm > 0.0005 Failed to execute (Select).

If i hard code that variable in a script as 0.0005 and run it as a stand alone it works fine.  Any ideas on why it won't run in arc?


texts:
I have thousands of tiff files with meta file containing corner coordinates. All images coordinates are different as they represent different locations.
How can I georeference all these files using corner coordinates in IDL or python?
Sample Meta file


texts: Symbol levels not working with Illustrator export

texts:
I use symbol levels to represent different overlaying lines by importance from a field value. The red lines are more important and yellow less. That's what I get in ArcMap:

And now that's what I get after an Illustrator export:

Symbol levels are not respected as yellow and orange lines overlay red lines.

I use ArcMap 10.2 (I don't know if the problem is still in 10.3).

I observed that this problem happens randomly. PDF and EPS exports are not concerned by this issue.

Here are the symbology and symbol levels properties for the layer:


texts: Keeping track of changes to MXD files using version control?

texts:
I am using ArcMap for the first time. My superior at the work place has asked me to do some research about the life cycle of an .mxd file.

Is there anything called a life cycle of an mxd file, like the life cycle of servlets?


texts: MapFish Print latitude, longitude innacuracy?

texts:
I'm new to GIS so this may be trivial: I'm using MapFish Print to create some PNG. I have issues to place the viewport. When I use a random GPS coordinate: -20.904522,55.475106 which seem to point roughly (<20m accuracy) to the same place in google and in openstreetmap.org, my viewport in MapFish Print is way up north (more than 100km).

What did I miss ?

Here's my config (I've commented some lines that I did try):

{
  layout: 'A5 landscape',
  title: '',
  units: 'dd',
  #srs: 'EPSG:900913',
  srs: 'EPSG:4326',
  #srs: 'EPSG:3857',
  #geodetic: true,
  #geodetic: false,
  outputFilename: 'mapfish-print',
  outputFormat: 'pdf',
  layers: [
    {
        type: 'OSM',
        format: 'image/png',
        baseURL: 'http://tile.openstreetmap.org/',
        ## This one seems incorrect and sets the 0,0 several
        ## thousand kilometer in the north
        #maxExtent: [-180,-90, 180, 90],
        maxExtent: [-180,-180, 180, 180],
        #srs: 'EPSG:900913',
        #geodetic: false,
        tileSize: [256, 256],
        resolutions : [
          1.40625, 0.703125, 0.3515625, 0.17578125, 0.087890625,
          0.0439453125, 0.02197265625, 0.010986328125, 0.0054931640625,
          0.00274658203125, 0.001373291015625, 0.0006866455078125,
          0.00034332275390625, 0.000171661376953125, 0.0000858306884765625,
          0.00004291534423828125],
        extension: 'png',
    },
  ],
  pages: [
    {
        center: [55.475106, -20.904522],
        #center: [0, 0],
        #srs: 'EPSG:4326',
        srs: 'EPSG:900913',
        #srs: 'EPSG:3857',
        geodetic: false,
        #geodetic: true,
        scale: 500000,
        dpi: 190,
        mapTitle: "",
        comment: "",
    }
  ]
}


texts: Get WMS version from QGIS layer in Python

texts: I need to retrieve the WMS version from a QGIS layer in Python. I can see the version in amongst the metadata(), but is there a better way than having to parse all of that HTML for one little value?


texts: I have a feature class that that I would like to dissolve. Upon the dissolve, I would like to concatenate a comment field, which is a text field. How can I do this? Conceptual idea below:


texts: Accessing portal for ArcGIS through web adaptor throws HTTP 404?

texts:
I'm in the process of configuring Web Adaptor with  Portal for ArcGIS . here's what happened . I Opened the Web Adaptor webpage then I was ordered to select a product to configure it with the Adaptor . I selected 'Portal for ArcGIS' , pressed next , then appeared a webpage with two lines :

Not Found

HTTP Error 404. The requested resource is not found.

This is exactly what I have done . The authorization process was successful and Portal for ArcGIS Windows service is running fine .


texts: Recording height profile from raster along a polyline shapefile?

texts:
I'm looking for the easiest way to record a series of raster values along a 400 meter transect in ArcGIS 10.3.1. I'm working with a LiDAR canopy height model with 3mx3m cells where each cell represents the height of the highest LiDAR return within that area. I'm hoping to produce an output similar to an elevation profile that shows the variation in canopy heights across a specific transect. I also want to show how the variations in canopy height are reflected by other layers generated through grey-level correlation matrix calculations (also a simple raster).

The goal is to create a chart that looks similar to this:

What is the easiest way to go about this in ArcGIS? I don't know anything about coding, but this seems like something that should be easy to do in Arc.


texts:
I have just started with Workflow Manager in ArcGIS Pro, and have created a Workflow Database, in which I can create new jobs.

However, I just see the three sample job templates (data edits, landbase updates and web edits), and was wondering how I can create my own job template / workflow.

Do I still need the old Workflow Manager desktop application, or can I do everything from within Pro?

I read the 'getting started' document but that's about running existing workflows (http://pro.arcgis.com/en/pro-app/help/workflow-manager/get-started-with-workflow-manager.htm)


texts: For loop within for loop with ArcPy cursor?

texts:
Below mentioned one is a part of my script.This is working when I assinged numbers for f.(Instead of "f" I have substitute 2.Then it is working).But i want to run this based on the no of rows of the table named "TL"(f).Can any one suggest a method to do this?

f=arcpy.GetCount_management(TL)
with arcpy.da.UpdateCursor(PCL_Info, ['PCL_PG_ID','PCL_IF_LN','PCL_IF_CN','PCL_IF_RM','PCL_IF_LU','PCL_IF_EL','globalid','created_user','created_date','last_edited_user','last_edited_date','adm_gn_ix']) as cursor:
    for row in cursor:
        for x in range(0,f):
           if(row[0] ==  parcel_ID[x]):
               c.append(row[0])
               c.append(row[1])
               c.append(row[2])


texts: Autofill raster bands when using batch processing in QGIS

texts:
If I have a single raster file with multiple bands that I want to do some processing on every band.
Is there a way of autofilling the bands when using batch processing in a similar way that you can add multiple files?
At the moment I have to manually add rows fill the filename down then select each band manually. This is fine for a layer with only a few bands but not if you have many bands.


texts:
I have made a mosaic of several old orthoimages. If it's possible, I'd like to do two things with the mosaic.

First item is to remove the edges of the smaller rasters that show up in the mosaic, so that the mosaic is more continuous looking. The edges show up as black strips and kind of break up the mosaic along where individual rasters lay on top of other rasters.

The second item is to try to even out the "color" in the images. These images are black and white and they come from scanning in the original b.&w. photos. Likely an artifact of the film, or maybe the scanning process (although I scanned all of these and did everything the same), one run of the orthophotos appears darker in the mosaic. I don't know if this is possible to fix. And if it is, if it won't mess up the rest of the detail of the mosaic.

While researching this question, another question pointed out tools in Arc that seemed useful. I don't know a lot about working with rasters (or GIS for that matter) and I tried these, but they don't seem to have changed much.

I am open on tools, Arc, R, QGIS, but I need some specific instructions if there is any Python code involved. Here is an image of the mosaic and an inset to point out what I mean about the uneven "color" of the mosaic.

I tried the methods suggested here and the raster images don't display. Instead, if I am at the full scale of the mosaic, 1:156,888, there is a framework, the "footprint", of the individual rasters, that kind of looks like a skeleton. If I zoom down to one raster, everything is the same shade of grey, no detail or features. FWIW, in the TOC, the Image has values between 2147483647 and -2147483648.


texts:
I have a raster which is a yearly time series, it is created with 30 raster. It is saved then as a raster with multiple bands.

I am trying to calculate the zonal statistics of all the bands at once. Imagine this like a "sandwich" (raster layers) and I want the mean and std of all the bands exactly where my finger press (zonal in vector file).

I used this analogy because, it is hard sometimes to express how the zonal statistics should be calculated for all the bands.

I have tried with ArcGIS, but I have the impression that the results are wrong. How can I do this either in QGIS or ArcGIS


texts: arcpy SearchCursor is creating a comma within a tuple

texts:
I am creating a list of tuples of coordinates of a feature class.

Feature class is created as such:

inFC = arcpy.GetParameterAsText(0)

The inFC is a point feature class.

I then have the following code to get the coordinates of all points in the feature class:

field = ['SHAPE@XY']
coordsList = []
with arcpy.da.SearchCursor(fc, fields) as cursor:
    for row in cursor:
        coordsList.append(row)

I was expecting an output of a list of tuples as below;

[(250000, 230000), (234567, 2859340)...]

However instead I am getting a list of tuples of a tuple and a commma, as below:

[((251075.0, 237048.0),), ((250687.0, 236717.0),),...]

What is the error in my code or my approach?


texts: Make OpenLayers request WMS tiles only within a certain bounds

texts:
Is there a way to actually tell OpenLayers not to request tiles out of particular bounds?

I want the rest of map and other layers still accessible, so maxExtent on map is not an option for me. I just want to restrict one WMS layer to particular bounds. If view port is inside the bounds and user zooms in fine, OpenLayers requests new tiles! but if viewport is outside, the WMS should only fill the bounded area with tiles from that WMS layer.

Does that make sense or is that even possible?


texts: how to keep relationship between feature classes after intersection one of feature classes?

texts:
I have geodatabase with 3 feature classes. one of them is circuit polyline feature class and two polygon feature class. substation and division

in geodatabase there is relationship between circuit and substation. I am going to intersect between division and substation. the result of relation  will be feature class with all substation and in attribute table all data about substations and division are available. but my problem is that there isnt the relation between circuit and substation in result feature class from intersect?
how can I make intersect feature class between sub and div that relation between substation and other classes maintain?


texts:
For a gone experiment I want to link the weather data to recorded track points. The 3 layers envolved are the track layer (points), the layer of privat weather stations ("PWS" points) and a geometry-less table of weather data (each row with "PWS_id" and "datetime").

I did already hard-coded join the track points with the most closest PWS by using the NNJoin plugin.

For each PWS I do have 1-n weather records per day.

Now I need to link/join the correct weather record to the track point by minimum "age" as difference to timestamp of trackpoint.

Is there a way to solve this problem from within QGIS (without python knowledge), or if postgres/-GIS is the best solution, what would be the SQL syntax to get the record with the minimum age(track-point, weather-record).


texts: Printing a multiline d.text in GRASS-GIS when using Python grass-session

texts:
Here is what I'm trying to do, but the text remain printed in one line. What is wrong?
from grass_session import Session
from grass.script import run_command

with Session(gisdb="path to mygisdb", location="myL72",
             create_opts="EPSG:31370") as sess:

    cptxt = "- first line \n"
    cptxt += "- second line \n"

    run_command('d.text',
                text=cptxt,
                color='grey',
                rotation=0,
                linespacing=1.8,
                at=[2, 45],
                align='ll',
                size=2.5
                )

But the result is - first line - second line.
Any idea?


texts: Align pixels from two rasters with different cell sizes

texts:
I have two rasters in a PostGIS database that I want to be able to compare. The first raster is at 5m pixel resolution and the second raster is at 30m pixel resolution and is on a different grid.

How I can align the rasters to the same grid based on the pixels, not on the rasters themselves?


texts:
Suppose I have a rather large NetCDF file that I open as a dataset using xarray.
import xarray

nc_in = r"....\D_P_WGS84_comp.nc"
nc_ds = xarray.open_dataset(nc_in)

It is composed of multiple raster variables, each with a single band, x, y coordinates and a crs (4326):
nc_ds

Suppose I want to visualize this very simply in my Python Notebook. I don't really want to add basemaps or legends or anything like that. I tried simply running show(nc_ds) but throws an errors AttributeError: 'Dataset' object has no attribute 'shape'`. I also noticed that the dimensions are frozen. What is the simplest way to visualize this dataset or a single variable of it in Python? Is there a mod or package extension that would do this?


texts: Calculating number of polygons or points that a set of polygons intersect with in ArcGIS

texts:
I have a set of ~100 polygons in ArcMap 10.2.2 that were generated using the minimum bounding geometry tool.

I have a second set of ~800 polygons generated the same way.  What I need to know is for each individual polygon in set 1, how many polygons from set 2 does it overlap with.

I know I can use a tool like 'intersect' but this will produce a list of all unique combinations of overlapping polygons which would be a hugely cumbersome output for a 100 by 800 data set.

Is there a way to get ArcGIS to simply return the number of polygons that each polygon overlaps with rather than a huge list of all unique intersections?


texts: Map put on hidden tab not shown after tab comes to front in OpenLayers 4?

texts:
I use OpenLayers4 with bootstrap-native.
My page has bootstrap tablist, and I put map layers on a tab which is initially hidden.
But after the tab is selected and comes to front, the map is white-out and is not shown..

How can I make it visible?

I tried

mapObject.renderSync();

but it doesn't work.

I also tried this map refresh but it throws error "TypeError: source.getParams is not a function".

I found a workaround by resizing the page, then map is shown.
But I want to display the map without having to resize the page.

How I can solve this?


texts:
I tried using the tool multiple times and my result values should be around 2000-3000 W/m2/d. However my results reach up to 6000 W/m2/d.

I am using a raster with cell size 30 meters and WGS 1984 UTM 48N. Sky Size is set to 400, Diffuse model type to STANDARD_OVERCAST_SKY and Diffuse proportion to 0,5 (since this is the actual diffuse proportion in this area ).

All other parameters in the tool are kept with their standard value.

Does anyone have ideas how to fix this issue?


texts: Building Interactive Map with Multiple Layers of Custom Shapefiles and Points data in Power BI

texts:
I want to built an interactive map that displays 2 types of data:

Polygons - I have 4 layers of custom shapefiles in ESRI
Points - lat/long data

Additional information:

I am limited to using Power BI for the front-end. Can use any other tools for the backend/GIS.
Needs to have drill-down capability into each layer.
Needs to be able to combine Reference Layers (Shapefiles) with Points data on one map.
Needs to include the ability to use multiple Custom Reference Layers (shapefiles)
Custom Reference Layers (shapefiles) can not be made Public (ArcGIS Online limitation)


texts: QGIS copy and paste error: "no features could be successfully pasted"

texts:
What does this error message mean, and how can you get more info on what might be the cause?

"no features could be successfully pasted"

(or more generally, not all of the selected items could be successfully pasted)

(I'm using QGIS 2.3 (build 422148c) on Vista Home Basic 32bit)

Seems like it's probably a problem with the source data; what's a good way to track that down?  If just modifying some of the source metadata would let the paste operation work, it would beat trying to reproduce the source data by hand in the target layer.

Thanks

UDPATE: looking at the metadata for successfully copied objects vs the object that dies, the one that won't copy has no firstX, firstY, lastX, lastY keys in the (Derived) section where the ones that paste successfully do have those keys with appropriate values.